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THERMAL WAVES IN IRRADIATED GRAPHITE 


By D. C. COOPER and N. E. HOSKIN 
(Atomic Weapons Research Establishment, Aldermaston, Berks.) 


[Received 9 October 1958] 


SUMMARY 
This paper presents the mathematical formulation of a problem in the release 
of Wigner energy in a graphite bar of variable specific heat and thermal! conductivity. 
The results of six such illustrative computations, obtained with the help of the 
I.B.M. 704 computer, are given to show the formation and motion of a thermal 
wave down the graphite bar. 


1. Introduction 


RECENTLY, Owing to the accident which occurred in the No. | pile at 
Windscale in October 1957, attention has been focused on the problem 
known as the release of Wigner energy. This problem arises because of 
an accumulation of energy in graphite which has been subjected to nuclear 
radiations over long periods (e.g. graphite used as a moderator in nuclear 
reactors). This energy (named after E. P. Wigner who first predicted the 


effect) is due to defects in the crystal lattice, and in nuclear reactors must 
be released periodically under controlled conditions, i.e. the graphite must 
be annealed. 

The Windscale incident, which occurred during such an annealing, re- 
sulted in the formation of a committee under Sir William Penney to 
investigate the cause of the accident, and a non-technical account of their 
findings has been published as part of a White Paper (1). Three other 
committees have also been set up under the chairmanship of Sir Alexander 
Fleck to study the wider aspects of the situation. The United Kingdom 
Atomic Energy Authority has been carrying out several investigations 
into the problem to assist these committees and the present paper is the 
result of some of this work. 

At the Geneva Conference held in 1955 several papers were presented 
under the general title ‘The Effect of Radiation on Reactor Materials’ (2), 
and recently a review of the effect of neutron irradiation on metals has 
been published by A. H. Cottrell (3) containing many references to the 
literature of the subject. For the purposes of this paper two papers from 
the Geneva Conference are relevant (4 and 5) as they contain information 
on the variation of the physical properties of graphite when it is irradiated. 
In the present paper we investigate the problem of the release of Wigner 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959) 
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energy in the form of a self-propagating thermal wave in order to determine 
whether such a wave is possible and if so to obtain information on the 
temperatures occurring and the speed with which the wave travels. 
Mathematically the equation describing the energy release is the partial 
differential equation of heat conduction with an extra term due to the 
generation of heat as the stored energy is released. If we consider the 
one-dimensional case of heat conduction along a rod this equation is 
c ( “4 oT a8 


of 2 SEO on 1 
— a OR (1) 


ex 
where the notation is explained below. The specific heat c and the ther- 
mal conductivity « are given functions of the temperature 7' and the 
amount of stored energy S, whilst the rate of energy release is assumed to 
be a known function of time and temperature. 

Analytical solutions of equation (1) have been obtained by A. J. E. 
Foreman if the specific heat and thermal conductivity are assumed constant 
and a simple expression is taken for the energy release term (6). These 
solutions show that a self-propagating thermal wave could form and travel 
down the rod, but usually a unique value for the speed of the wave is 
not obtained. It is the purpose of this paper to describe a method of 
solution of equation (1) which has been programmed for computation by 
the I.B.M. 704 computer, the programme being written so that any of 
the parameters of the problem may be adjusted and solutions quickly 
obtained of problems whose basic equations are the same. The specific 
heat and thermal conductivity of the bar may be functions of the tem- 
perature, or of the stored energy, or of both these quantities, and the 
programme is not restricted to a constant initial temperature and stored 
energy distribution. The machine time required is about 8 minutes for 
every hour of real time computed. 


Notation 
xr distance along the bar from the hot end (cm). 
t time in sec. 
T Temperature. 
T,, and T, temperatures of the sources at the hot and cold ends. 
T, initial temperature of the bar. 
hy, and h, heat transfer coefficients at the hot and cold ends in cals 
°C-1 cm-* sec", 
specific heat (at constant pressure) in calories per gramme 
per °C. 
thermal conductivity in cals °C-! em-! sec-!. 








THERMAL WAVES IN IRRADIATED GRAPHITE 395 


frequency parameter such that S,d(Inv) is the amount of 
energy stored in the range d(Inv) at Inv. 
total stored energy in calories per gramme. 
density in grammes per c.c. 
activation energy in eV. 
Boltzmann’s constant. 
B = Q/k. 
All other symbols used are defined as they arise. 
All temperatures are measured in degrees absolute unless otherwise 
stated. 


2. Energy release kinetics 

At the present time the theory of the method of release of energy in 
irradiated graphite is incompletely understood. It has however been 
pointed out by W. M. Lomer that there is a simple model which may be 
useful for predicting the rate of energy release along an arbitrary tem- 
perature-time curve which is not very far from that found in some previous 
experiment, in which the rate of energy release was measured. This is 
what may be called the constant activation-energy model and the follow- 
ing description by Lomer is to be found in a modified form in (7). 

We assume that the stored energy is released on annealing as if it were 
a range of first-order processes, all with the same activation energy. These 
processes may be characterized by the frequency v (or more conveniently 
Inv) which appears in the rate equation 


dS, 
dt 


where Q is the activation energy and k is Boltzmann’s constant. The total 
stored energy per unit mass is the sum of all these contributions 


= —vS,exp(—Q/kT), 


s= f Sdn»), 


—-@® 


Integrating the first of these equations we have 


t 
S,= Syexp| — ( vexp(—Q/kT) at 
és 


where S, is the initial density of the energy stored at the frequency rv, this 
initial time being referred to as t,. These equations provide the necessary 
information to calculate the rate of energy release if the initial stored 
energy spectrum S, can be obtained from some previous experiment, but 
may be simplified because of the sharp dependence of the exponential 
function on temperature. The quantity S, passes very rapidly from S, to 
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zero, this sharp disappearance occurring at about a frequency v, given by 


t 
V1 ( exp(—Q/kT) dt = 1 
i 
and so the total stored energy may be written as 
Iny, 
S= ( S,d(Inv). 

Dropping the suffix 1 we therefore have the following two equations to 
determine the stored energy function in terms of the temperature de- 
pendence at a given point: 

Inv 
S = | S,d(Inv), (2) 
where S, is a given function of Inv and 
t 
= | exp(—B/T) dt. (3) 
ti 
The constant f, corresponding to Q/k above, has been taken to be given 
by f = 17,980 degrees absolute, 
which corresponds to an activation energy (Q) of 1-55 eV. 

It must be noted that in equations (2) and (3) the temperature has to 
be expressed in degrees Kelvin. The initial stored energy spectrum used 
(S,) has been obtained for a particular specimen from a laboratory experi- 
ment at Harwell and is tabulated in Table 1 together with the correspond- 
ing value of S determined from equation (2), which is also plotted in Fig. 2. 
It is used here for illustration and is not necessarily the best fit to all the 
experiments. 

The initial value for v is obtained by making the assumption that the 
graphite has been brought up to its initial temperature 7, at a given 
constant rate (d7'/dt) and can be obtained as follows: 

Write u = 8/T in equation (3) and use the fact that d7'/dt is constant; 
then 





OBR? ioe { exp(—u) 5, _B (| ~ (4) 
v dT’ /dt u* dT /dt u® a 
where, since u is of order 40, we have used the asymptotic formul. 
(“e pee exp(—z) 
u* es: aps 


z 
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TABLE 1 





Se S+constant In v Y S+-constant 





0°54 o6 30 38°63 
o-72 1°26 ° 41°1275 
0°93 2-91 31 3: 43°745 
118 5°02 “s 40°4575 
1°33 6-275 32 ; 49°2475 
1°47 7°675 32°5 . §2°0725 
1°73 9°275 33 ; 54°8875 
1°94 10°1925 33°5 . 57°6475 
2°05 Irtg 34 5 60°2875 
2°07 12°22 35 5 64°9775 
2°06 13°2525 36 P 68-6075 
2°06 14°2825 37 “ 71°0875 
2°09 15°32 38 72°6975 
2°19 17°46 39 73°6475 
2°31 19°71 40 749175 
2°46 22°095 41 74°3425 
2°64 24°645 42 744425 
2°90 27°415 43 , 74°4975 
3°30 30°515 44 , 74°4975 
4°04 347185 50 ’ 74°4675 
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By ‘final’ in equation (4) we refer to the state to which the graphite has 
been brought before we start our energy-release calculations (i.e. the tem- 
perature 7). Because of the rapid variation of the exponential function 
we can neglect the value of exp(—w) in whatever configuration is taken 
as initial and thus obtain 


T 
v0 = FE PPT) (5) 


where v, is the starting value to be taken for v corresponding to the initial 
temperature 7}. 


In our computations we have taken 


qr = 2 © per min 


which gives as our starting value for v 
(i) 7, = 123° Inv, = 39-83, 
(ii) 7, = 100° Inv, = 42°76. 


3. The boundary conditions 

The physical model used is shown in Fig. 1. The problem is treated one 
dimensionally, a graphite bar of length 1 metre having heat sources S,, 
and S, placed close to either end whose temperatures can either be kept 
constant or be given functions of the time. The wave is generated by 
raising one of the sources (S,,) to a temperature above that of the other 
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(S,). Heat is assumed to transfer from the sources to the bar (or vice- 


versa) by means of the Newton law with constant coefficient, i.e. we have 
eT 
hy(T,—T) = —«— at the hot end | 
éa 


(6) 
and h(T.—T) = +0 at the cold end 
a 
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where the notation is explained in section 1. The distance variable z is taken 
to increase on going from the hot to the cold end of the bar. 

In the actual examples taken the source S,, has been made to increase 
from the initial temperature of the bar (7)) to its final temperature (7;) 
according to the law 


Ty = Ty+(T,—Ty)exp(—t?/A*), (7) 


where A is a constant which is chosen so that the final temperature (7;) is 


achieved in about half an hour. The values taken in the solutions reported 
here are 


A 10 (min-'), 7, = 400°C, TT, = 100° C and 123°C. 


The temperature of the source at the cold end (7) has been kept constant 
and equal to the initial temperature of the bar and any heat loss from 
the sides of the rod has been neglected. 


4. The physical variables 


The formula for the variation of specific heat with temperature given 
by Magnus, and quoted in (8) has been used. This is 


C,, = 7-635 4 39-06 x 10-87’ — 43-02 x 10-672 + 


+ 29-57 x 10-®73— 11-01 x 10-74, (8) 


where C,, is the specific heat at constant pressure measured in joules per 
gramme atom per °C and T is the temperature in degrees centigrade. 
To obtain the units for this problem we use the conversion formula 


1 
a ae ee 
. 50-16 ”” 


where c is the specific heat in calories per gramme per °C. 
The thermal conductivity is a function of the stored energy whose 
variation can be represented by the formula 
Ko 


~ 4S+ constant 
K 


with “o _ 19-5 initially 


K 





and Ky = 0°33 cal °C-! cm-" sec! / 


where S is the total stored energy measured in calories per gramme 
(see, for example, Kinchin (5)). 
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The density of the graphite used is given by 
p = 1-65 g per c.c. 


and the Newton coefficient for heat transfer at both ends of the rod is taken 
to be 


— J cal °C-! em-? 1 

hy = h. = gy cal °C-' cm-? sec. 
Since there is some uncertainty about the value to be used for the last 
constant, the calculations have been performed three times, once with 


the above value, once with four times this value, and once with one-quarter 
of this value. 


5. The finite difference representations 

(a) The heat conduction equation 

The finite difference representation of equation (1) is obtained by 
considering the bar to be split up into NV equal sections each of width Az, 
at any time the temperature and stored energy being given at the NV +1 
ends of each section and also at the N midpoints of each section. In each 
section we assume a parabolic distribution of the temperature and stored 
energy, the parabola being uniquely determined by the given values at 
the end and midpoints of the section. Equation (1) is then satisfied in an 
integrated form through each section, and conditions are also imposed to 
ensure that the heat flow out of one section is equal to the heat flow into 
the next section. 

The notation to be used for the finite difference representation is that 
subscripts 2n, 2n+-1, and 2n+2 are used to indicate the left-hand end, 
the midpoint, and the right-hand end of the nth section whilst super- 
scripts ¢ and ¢+-1 are used to denote two successive time steps. 

Integrating equation (1) with respect to x through the nth section we 
have 


2n+2 2n+2 
ft a oe ee é 
K - p- éT dz = p— S dx, 
ox |, et . ot 
_ 2n Qn 


where ¢ is used for the mean value of c through the time interval being 
considered. 

Using the assumed parabolic distribution for S and 7 through the 
section, and replacing time derivatives by differences, we obtain the finite 
difference representation of the heat conduction equation for the nth 
section in the form 


A Tr + BTS 1 t CTtt1 + fc 2 | 4St+1 


2n+2 + Soi} 


2n+17 2n+2 


DT n+ ET in +17 FT 4n42+ Stn+45in +17 Stn+2 (10) 
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A = €y,—3ax'gn} 2—9axyy 
“= 4¢,,, +1 -t- 12axft? + 12ax4t! 
= Con _— Panis -—3anhi? 
; €a,+ Sakon 42+ 9ar'y 
= 4sn41- 12aktn +27 12ar'yn 
Conse t Daron 4 9+ Sauce, 
1 At 
a= -——>: 
p (Ax) 
There is also the condition for continuity of heat flow between sections, 
i.e. the quantity at 
Ox 
must have the same value at the right end of the nth section and the left 
end of the (n+-1)th section. 
In finite difference form this condition is 
37 2, +27 4T,,, +1 + T 2» T 2», 14473, +3 37, +2° (1 1) 


Equations (10) and (11) provide an implicit system of equations which 
may be solved by an iterative method in order to progress the computa- 
tion from knowing all variables at one time to knowing all variables at 


the next time step, if we also have equations to find the energy release S. 
The system is stable for any values of the space and time intervals. 


(b) The energy release equations 

To find the cut-off frequency v equation (3) has to be put into finite 
difference form, but care must be taken because of the rapid variation 
of the exponential function, 

We replace equation (3) by its equivalent 


d(iInv) = —exp(Inv—f/T') dt, (12) 


and, in progressing from one time step to another, take the mean of the 
expression in brackets on the right-hand side of equation (12). This avoids 
computing with very large and rapidly varying numbers since both Inv 
and 8/T are of the order 40. 

‘quation (12) enables a new value of Inv to be found but we cannot 
use it directly since In vy occurs on the right-hand side. If on the right-hand 
side of (12) we use the previous iterate to Inv we obtain d(in v) and hence 
a new Inv; we have thus an iterative scheme which could be incorporated 
in the iterative scheme for finding 7 in the last section. This process has 
turned out to be divergent in certain cases (in particular if we keep T 
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constant in equation (12) the process can be proved to be divergent if 
A(inv) > 2). 


Instead of using equation (12) in this manner we therefore use 


A(inv) = > (A In») —2}, (13) 
where E = Atexp mean(Inv—8/T), 


and an accent refers to the value of the quantity at the previous iterate. 
This formula is obtained by applying the Newton—Raphson process to 
find the root of equation (12) regarded as a transcendental equation in 
Inv, neglecting variations in 7’, and gives a convergent process even for 
large values of Alnyv. 

Having obtained Inv we can find the total energy release S at the point 
from equation (2). The quantity S, is a given function of Inv and we can 
therefore perform the integration of equation (2) independently from the 
rest of the calculation and provide a table of values of S corresponding 
to Inv values, apart from an irrelevant arbitrary constant in S (see Table 1). 


(c) The boundary conditions 


Using the assumed parabolic distribution of temperature in the first and 
last sections of the rod we can replace equations (6) by their finite difference 
representations as follows: 


(34 °)n, = “a 7,+ 47,7, 


Ko Ko (14) 
(hh. Aa h, Ax 
3-4. ~\ ex = - i+ $Tyy-1—Tav-a | 
\ Kon Kan 


where 7;, and 7, are the temperatures of the two heat sources and are 
given functions of time. 


(d) Interval size 


The intervals used for all the runs were 
At = 20sec and Az = 5cm. 


The time interval in one of the cases (initial temperature = 123° C and 
h = 1/900) has been halved, with no resultant change in the speed of the 
wave and less than | cm change in the wave position at any time. 

The space interval in the same case was also halved, the results of this 
being shown in Fig. 4. It can be seen that the irregularities near the wave 
extremities have been evened out and that the wave slope is slightly 
greater, but there is no change in the wave speed. 
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—— initial temperature of 123°C 
i) 
—<<= Initial temperature of 100 C 
Heat transfer coefficent in calories 
°C em" sec™ 
h=-00111 Ah=-00444 
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€ 
S 
) 
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73 
O 





Time (hrs) 


Fic. 3 


6. The steps in the computation 


A summary of the steps necessary to advance the calculations for each 
time interval is given in this section. We assume we know all quantities 


throughout the rod at some time ¢ and wish to find the quantities at time 
t+-At. 


(1) Obtain approximations to the temperature and Inv at all points of 
the rod at time t+ At. 

(2) Use equation (13) at all points of the rod to find better values for 
Inv at time ¢+-Af. 

(3) Find S from these values of Inv by reference to a previously con- 
structed table of S as a function of Inv. 

(4) Go through the rod using equation (10) to predict better values of 
T at the midpoints of each section and equation (11) to predict 
better values of T at the end points of each section. 
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(5) Use equations (14) to determine better end temperatures of the rod 
by means of the boundary conditions. 

(6) Compare all these temperatures with those of the previous iteration 
and repeat from stage (2) if the process has not yet converged to 
the desired accuracy. 
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A test is inserted so that points of the rod which have already converged 
and which cannot be altered by any unconverged points are not recom- 
puted at every iteration. 


7. The results of the computations 

As previously stated, illustrative calculations have been performed for 
two different constant initial temperatures of the bar, each of these 
calculations being carried out for three different values of the coefficient 
of heat transfer. The results are summarized in Table 2, the formation 


TABLE 2 





Initial 
temperature 


Time taken 
to form wave 


Final wave 


speed 





(min) 


(cm/min) 


00002778 32 21 
O-OO1IIII 17 2°0 
00044444 10 2-0 


0:0002778 7 ee 
O°OOI iit 21 12 
00044444 Ir 1°2 














and motion of two of the waves is plotted in Figs. 4 and 5 and a graph 
of the position of the wave in all six cases is plotted against time in Fig. 3. 

The results show that in all these cases a thermal wave was formed, 
which then travelled down the bar with a constant speed dependent on 
the initial temperature of the bar but independent of the heat transfer 
coefficient, i.e. of the way in which the wave was formed. The time taken 
for the wave to form, however, increases as the heat transfer coefficient is 
decreased. 


8. Acknowledgements 

Throughout the whole project much advice and assistance has been 
received from Mr. E. P. Hicks. Dr. W. M. Lomer of Harwell provided 
help with the physical formulation of the problem. 


REFERENCES 


1. The Accident at Windscale No. 1 Pile on 10th October, 1957, H.M.S.O. Cmd. 302. 

2. Proceedings of the International Conference on the Peaceful Uses of Atomic Energy 
(Geneva 1955), vol. 7, section 11s. 

3. A. H. Corrrett, ‘The effect of neutron irradiation on metals and alloys’, J. British 
Nuclear Energy Conference, 3 (1958), no. 1, p. 50. 








406 THERMAL WAVES IN IRRADIATED GRAPHITE 


4. 


DIMM 


W. K. Woops et al., ‘Irradiation damage to artificial graphite’, see reference 2, 
455. 

G. H. Kincutn, ‘The effects of irradiation on graphite’, see reference 2, 472. 

A. J. E. Foreman, A.E.R.E. Theoretical Physics Memo. T/M 162 (1958). 


. W. M. Lomer and D. E. Rimmer, A.E.R.E. Report M/R 2603 (1958). 
. International Critical Tables, vol. 5, p. 94 (1929). 





DIFFUSION FROM AN INSTANTANEOUS POINT 
SOURCE WITH A CONCENTRATION-DEPENDENT 
COEFFICIENT 


By R. E. PATTLE 
(Ministry of Supply, C.D.E.E., Porton Down, nr. Salisbury, Wilts.) 
[Received 24 July 1958] 


SUMMARY 
An expression is given for the concentration distribution produced by diffusion 
from an instantaneous point source in one, two, or three dimensions, when the 
diffusion coefficient varies as a positive power of the concentration. 


1. Introduction 

THE equations of diffusion, and the corresponding equations of heat 
conduction, occur in a number of physical problems; in some of these, 
however, the coefficient varies with the concentration of the diffusing sub- 
stance. There is reason to suppose, for instance, that the seepage of liquid 
into a porous body or cloth, under the influence of capillary forces, may, 
under certain conditions, be described in terms of the diffusion equation; 
in this case, the seepage coefficient (the analogue of the diffusion coefficient) 
increases rapidly as the amount of liquid in the pores increases, and the 
resistance to viscous flow decreases. 

The present paper describes a particular case of such diffusion. 


2. Diffusion from an instantaneous point source 

If a quantity Q of a substance (quantity being defined as the line, area, 
or volume integral of concentration, as appropriate), is liberated at time 
t = 0 at a point in an infinite line, plane, or volume, and diffuses with a 
constant diffusion coefficient D, the concentration at any point at distai:ce 
r from the point of liberation, at time t, is given by (1) 


C = Q(4rDt)-¥ exp(—r*/4Dt), (1) 


where s is the number of dimensions (1, 2, or 3) in which diffusion is taking 
place. 

A solution of the equations of diffusion is now presented for the case 
of diffusion of a quantity Q from an instantaneous point source when the 
diffusion coefficient varies as a positive power of the concentration, being 


aww D = DiC/C,)" (n> 0). (2) 
The equation for outward diffusion from a point in s dimensions into 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959) 
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a medium of zero mean motion may be written 
eC/et = r°-%e/ar{r*-» D(eC /ér)}. (3) 
A solution of this equation appropriate to an instantaneous source is 
C = C,(t/ts) 8" {1 —r2/r2}!™ (where t > 0, and r? < rj) = (4a) 
and C=0 (where r? > r?). (4b) 


Here the symbol r, (the radius of the zone of diffusing substance), is 
A 1 
written for convenience for the expression 


+2) 
r{t/t,)o"t® = ¢,. 


The quantities ¢, and r, occurring in equations (4a) and (4b) are inde- 
pendent of distance and time and linked to the originally postulated 
quantities Q, D,, and C, by the equations 


r§ = (Q/C,)nr-*{T(1/n+ 1+ $8)/T(1/n+1)} (5) 
and ty = r2n/2D,(sn+ 2). (6) 


If C is suitably differentiated with respect to r and ¢, it can be shown 
that equations (4a) and (4b) are solutions of equation (3), subject to the 
condition of equation (6). Equation (5) is derived from the constant 
quantity condition. In bringing the one, two, and three-dimensional cases 
into a single formula, the identity 

(1, 2, 27) wi#/T'(48s) (8 = 1, 2, or 3) 
may be used. As n> 0, equations (4a) and (4b) tend to the form of 
equation (1), which is valid when D is constant. 

The equations (4a) and (4b) represent a region of diffusing substance 
with a definite edge, specified by the boundary between (4a) and (4b); 
the diameter of this region (2r,) varies as the 1/(sn-+-2) power of the time, 
and the concentration at the origin varies inversely as the s/(sn+-2) power 
of the time. For instance, if the diffusion coefficient is proportional to 
concentration, and diffusion is taking place in a three-dimensional space, 
the substance is contained in a spherical volume whose diameter varies 
as the fifth root of the time elapsed since the liberation of the substance; 
ithe concentration at the centre varies inversely as the 0-6 power of the 
time. If the seepage coefficient in a sheet of blotting paper is proportional 
to the square of the concentration of ink in the paper, and a fixed quantity 
of ink is liberated in a sufficiently small region of it, the diameter of the 
resulting blot will, according to equation (4a), vary as the sixth root of 
the time elapsed since the liberation of the ink. 

When n > 1, the concentration gradient at the boundary of the region 
is infinite; when n = 1, it is finite, and when n < 1 it is zero. Whenn = 1, 
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the graph of concentration against distance along a line through the origin 
is a parabola; when n = 2, it is one-half of an ellipse. In the example 
given above, the blot will have a sharp edge. 


+0 





i i 


05 0 
r/\CofQ) 
Fic. 1. Concentration C due to diffusion in two dimensions from ap 
instantaneous point source of quantity Q of a substance whose diffusion 
coeflicient varies as the nth power of the concentration. Abscissa: distance 
from origin in units of .(Q/Cy). Ordinate: ratio of concentration C to 
concentration (Cy) at centre. Curve 1, n = 2; Curve 2, n = 0-5; Curve 3, 
D constant. 





1 
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Fig. | shows the distribution of equal quantities Q of a substance, 
diffusing in two dimensions, at the instant when the concentration at the 
origin has fallen to Cy. Two of the curves represent the cases where the 
diffusion coefficient varies as the square and square root of the concentra- 


tion, and the curve of (1) (for constant diffusion coefficient) is given for 
comparison. 
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SUMMARY 

The well-known solution for laminar forced convection in a uniformly heated pipe 
neglects the effect of buoyancy forces caused by temperature variations in the 
fluid. This effect depends on the orientation, but in a horizontal pipe it produces 
a circulation of the fluid in a direction normal to the pipe axis, with a consequent 
modification of the main flow. The present treatment includes these buoyancy 
effects, but is restricted to small rates of heating (which correspond with small 
temperature gradients along the pipe wall) so that the motion due to buoyancy can 
be regarded as a secondary flow. In this way solutions for the velocity and tem- 
perature fields far from the pipe entrance have been obtained as power series, which 
are shown to depend on the dimensionless parameter AB, where A is a Rayleigh 
number based on temperature gradient along the pipe wall and B is a Reynolds 
number based on pressure gradient along the pipe axis. The Nusselt number and 
the resistance coefficient have been calculated up to the terms in (A B)*, and the 


forced convection solution is shown to be in error by about 10 per cent for 
AB 3,000, 


1. Introduction 

Tue forced convection solution for laminar flow along a straight pipe 
with walls maintained at a uniform temperature gradient (which corre- 
sponds with uniform heating) was given by Nusselt, and may be found 
in Goldstein (1, 622). This solution is entirely independent of the pipe 
orientation since gravitational forces are ignored; it predicts a Nusselt 
number of 6 based on pipe diameter and mean temperature across a 
section, and the same volume flow and wall stress as for normal Priseuille 
pipe flow. 

In this type of forced convection solution the fluid motion is regarded 
as due entirely to external forces (in this case the pressure gradient forces) 
and the heat transferred is assumed to diffuse through the fluid and be 
swept along by it without affecting the local density in any way. Thus 
forced convection represents a limiting approach for vanishingly small 
heat transfer, and these solutions are valid approximately only within a 
very limited range of flow and heat transfer conditions. Moreover the 
solutions do not, in themselves, permit any estimate of the extent of this 

[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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range within which they are physically realistic, and previous discussions 
of the effectiveness of solutions have of necessity been based on experi- 
mental observations and not on a comparison of characteristic values of 
the inertial and buoyancy forces for the particular system. 

One way of investigating probiems of convection in which the applied 
and buoyancy forces are comparable is to use the forced convection 
solution as a first approximation to both the velocity and the temperature 
fields. The next approximation can then be found as a secondary flow 
driven by the buoyancy forces and a secondary temperature distribution 
due to the modified field of flow. This merely gives a perturbation of the 
forced convection solution, but it is sufficient to yield some information 
on mixed convection; and it does allow for interaction of the velocity and 
temperature fields, which is the essential feature of thermal convection. 

This procedure will now be applied to laminar convection in a long 
horizontal pipe, the walls of which are heated at a uniform rate. Solutions 
will be obtained as power series in the Rayleigh number A, and will be 
convergent for sufficiently small values of A. The secondary flow in a 
curved pipe (treated by Dean 2, 3) and that in a pipe rotating about an 
axis perpendicular to its own (treated by Barua (4)) are broadly similar 
to the flow treated below, and similar methods of analysis will be used. 


2. Equations for convection in horizontal pipes 


Consider the steady laminar motion of a fluid along a horizontal circular 
pipe of radius a, the walis of which are heated uniformly so that a constant 
temperature gradient 7 is maintained in the direction of the axis. The 
flow will be referred to cylindrical polar coordinates (r,¢,2), with ¢ 
measured from the upward vertical and the z-axis along the axis of the 
pipe; the velocity components are denoted by u, v, w. It will be assumed 
that the effects of dissipation and of the pressure term in the energy 
equation can be neglected, that variations in density due to temperature 
changes are so small that they need be taken into account only in the 
buoyancy term, and that the kinematic viscosity, v, and the thermometric 
conductivity, «, can be taken as constants. Since the wall temperature 
increases slowly with distance these assumptions are reasonable when 
applied to regions of limited length in the pipe; they will introduce some 
quantitative errors into the solution, but they should not change its 
general character. 

The equation of continuity is then 

oru ev Orw 
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the energy equation is 


t z oa Lw or = KV2T, (2) 


or =r ad ox 


28 ‘ os 
F183. 2b?» # 
where = —34+--=-+5mmte 
Or? rodr’ r® Gd? Ox? 
and T7' is the local fluid temperature. The momentum equations can be 
written 


‘ a 9 a > De 
ou vou ou si 1 Op ; u 2 oO ‘ 
6 son he as ee he Pe Py (02 -——-—- Bg(T,,,— T )eos 4, 
or r ch Cx r peor \ ie r? od 
(31) 
ov. vw Ov Cv uw l Op ? 20u ov ioe 
6 —+.— — +o —-+4 ee o +vi V%v+— — — ) +- Bg(T,,.— T)sin ¢, 
or red ae pr od rad r* 
(3ii) 
Cw . v Ow ow l op — 
4 —+-—+w. —_- Pi ww, (3 iii) 
or  r Op Ox p 0x 


where p is the local fluid density, 8 the coefficient of expansion for the 
fluid, g the acceleration due to gravity, and 7,, the local temperature of 
the pipe wall. In equations (3i) and (3ii) the buoyancy force has been 
calculated relative to fluid at the same level adjacent to the pipe wall and 
the remaining distribution of force has been absorbed into the pressure, p. 

For steady convection sufficiently far from the pipe opening to avoid 
inlet length effects, the temperature throughout the pipe increases uni- 
formly with distance along the axis. Hence the distribution of buoyancy 
force in sections of the pipe is independent of x; and, as the secondary 
flow is caused by the buoyancy forces, the flow field must also be inde- 
pendent of x. It follows that there should be a similarity solution with 
u, v, w, and T,,—T as functions of r and ¢ only, and with p = ya+ P(r, ¢) 
(where P contains the terms absorbed into p). 

The continuity equation reduces to the form 


hence a dimensionless Stokes stream function % can be introduced in such 
a way that 
° ruly = eb/ad, v/v = —op/or. 
The main equations (1-3) may be reduced to non-dimensional form by 
the transformations r = aR, w = (v/a)W, and 7,,—T = rao, where o is 
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the Prandtl number (v/«x). If the pressure is eliminated between equations 
(3i) and (3ii) the resulting momentum equations are 


141% © Flory a(S rings} Mooag), 
Viet Ae ob aa) td flee aR in t+ Bagot “) 


ews l(ab @ eb O\y yp 
ve W+ Gea 34 ap) +48 = 0. (5) 


ae 
5 ait at a ya" 
where A = Sgra*/xv is the Rayleigh number, B = —ya*/4pv? = ore 
is the normal Reynolds number for Poiseuille pipe flow based on the pipe 
diameter and the mean velocity across a pipe-section, and 
». & Beet # 
1 oR ROR” Ri ag™ 
The form of the energy equation depends on the temperature boundary 
condition at the wall of the pipe. Pipes will usually have reasonably thick 
walls of material with thermal conductivity much higher than that of the 
fluid, and so there will be little variation in wall temperature round the 
pipe girth, and this will be specially so for slow rates of heating when 
asymmetries of the flow will have small amplitude. Hence write 7, = T,+72, 
where 7), is the wall temperature in the section containing the origin. 
Although there is uniform heat transfer per unit length of pipe, the local 
heat transfer will be slightly greater near the bottom of the pipe than near 
the top. Equation (2) now reduces to the non-dimensional form 


+W=0. (6) 


The boundary conditions are: 
u, v, w, and @ are zero on r = a; 
u, v, w, and @ are finite for r = 0. (7) 
Although solution of equations (4-6) satisfying conditions (7) is a 
matter of considerable difficulty, successive approximations to the solution, 


can be determined by expanding ¥, W, and @ as power series in the Ray- 
leigh number, provided that this is numerically small. Accordingly 


suppose that b == Ay, + A%e+... 
W = W4+AW+A*W+..., 
= 0,+A0,+AO,+.... (8) 


The leading term yf, of 4 must vanish because there is no circulation when 
A is zero, but 4, is not zero since the difference in temperature between 
a fluid element and the neighbouring wall is proportional to 76. When 








414 B. R. MORTON 


relations (8) are substituted in equations (4), (5), and (6) three sets of 
equations for the functions ,;, W,, and 6; are obtained by equating coeffi- 
cients of powers of A. 
From (5) the basic equation is 
V?W,+4B = 0, (9) 
and the solution W, B(i — R?) (10) 


satisfies the boundary conditions W, = 0 at R = 1, W, finite at R = 0. 
This is ordinary Poiseuille flow under a pressure gradient —4pv?B/a? in 
an unheated pipe. 

The first of the equations derived from (6) for the temperature distri- 
bution is v26,+W, = 0 (11) 
subject to 6, = Oat R = 1, 4, finite at R = 0, and the solution is axially 
symmetrical as no account has been taken of gravity at this stage of the 
approximation. This is satisfied by 


0, = 7 B(1— R*)(3— R?), (12) 


which is the customary solution for forced convection. 
The first order approximation for y satisfies the equation 


Vid, Crosin (13) 


obtained as the coefficient of A from equation (4), and the boundary 
conditions éy,/8R, eb,/6d = 0 on R = 1, and R-"(éb,/6d), &b,/eR remain 
finite at R= 0. If y, is assumed to have the form »%,(R)sind, the de- 
pendence of equation (13) on ¢ is eliminated, and y,(R) is easily found, 
yp, jpn BD R(1 — R*)?(10— R?)sin ¢. (14) 
Similarly, the first-order approximation to W satisfies the equation 
1 ab, eW, 
R 0d @R’ 





ViW, = (15) 
and the boundary conditions, W, = 0 at R = 1 and W finite at R = 0. The 
dependence of equation (15) on ¢ is eliminated by taking W, = W,(R)cos¢, 
whence 


W, sajysp B? R(1 — R*)(49—51 R2+ 19.R*— R®)cos 4. (16) 
The first-order approximation to @ satisfies the equation 
a 00 
VP, = 2 ey 17 
7 se (17) 


and the boundary conditions, 0, = 0 at R = 1 and @, fiite at R = 0. 
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Hence 
0 


| = —sortsaoo B2R(— R2){(381 + 13250) — (354+ 16750) R?+ 
+ (146+ 9250) R*— (29+ 2000) R®+ (1+ 20) R8}cos¢. (18) 











r 


5 


7 





L1-0 
Fic. 1. Fic. 2. 


Fic. 1. The streamlines corresponding to the first-order approximation for 
flow in planes normal to the pipe axis. (Only the right-hand half is shown.) 


Fic. 2. The first-order approximation to the velocity components in planes 

norma! to the pipe axis plotted in non-dimensional form for the case AB = 1,000. 

The radia! component U is a profile taken vertically through the pipe axis, and 

the transverse component V is a horizontal! profile through the axis. The second- 
order approximation does not affect these profiles. 


This completes the approximate solution to the first order in A for 
convection in a horizontal pipe which is heated so that there is a constant 
temperature gradient along the walls and uniform temperature round the 
girth. The streamlines corresponding to this secondary flow in the right- 
hand half section of the pipe are shown in Fig. 1. Fig. 2 shows the com- 
ponents of the non-dimensional velocity in a section of the pipe normal 
to the axis for the particular case AB = 1,000 and o = 0-73 (for air); 
the radial component U is plotted for a vertical traverse through the pipe 
axis (i.e. 6 = 0 for the top half and ¢ = a for the bottom half), and the 
transverse component V for a horizontal traverse through the axis (i.e. 
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¢ = 4m and 37/2); in each case the first-order approximation is shown by 
the continuous line (note that U = au/v and V = av/v). Fig. 3 shows the 
longitudinal velocity plotted as W/B and the temperature difference 
plotted as @/B (with zero displaced to the right to avoid overlapping of 
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Fic. 3. The fluid velocity parallel to the pipe axis and the 

temperature difference plotted non-dimensionally as W/B and 

0/B, respectively, for AB 1,000 and o = 0-73. The R-scales 

have been separated for clarity. The dotted line is the forced 

convection solution, the continuous line includes the first- 

order corrections and the dashed line includes also the second- 
order corrections, 


the curves), in each case for a vertical traverse through the axis and with 
AB = 1,000 and o = 0-73; in this figure the dotted lines show the zero- 
order approximation of forced convection, the continuous line includes the 
first-order correction and the dashed line includes the second-order as well. 
These diagrams are of interest because they give some indication of the 
way in which the pattern of flow will develop with further increase of the 
parameter A B. 

It may be noted from equations (8, 10, 12, 14, 16, 18) that the full 
convection solution depends essentially on the product A B of the Rayleigh 
and Reynolds numbers.t The reason for this is clear if it is recalled that 


+ This parameter AB can be regarded as the ratio (characteristic inertia force)*/(charac- 
teristic buoyancy force)(characteristic viscous force). 
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A is proportional to the increase in wall temperature per unit distance 
along the pipe and B is proportional to a characteristic velocity of flow 
along the pipe. An increase in A means that the fluid will be carried 
through larger temperature differences and hence there will be increased 
buoyancy forces; but the same effect can be produced by an increase in 
B, and hence in the velocity of the main flow. 

The two most interesting consequences of such a solution are the modi- 
fication of the heat transfer due to the secondary flow, and the effect of 
heating on the tangential stress at the wall (or on the rate of flow for given 
pressure gradient). The latter effect arises because the secondary flow 
carries fluid from the centre of the pipe, where the velocity is high, towards 
the walls, where the velocity is reduced by the wall stress, with consequent 
loss of energy. However, as both W, and @, vary with cos¢, neither the 
flux of fluid nor the heat transfer across the whole section of the pipe will 
be changed by the first-order approximation, although there will be local 
variations. 

To find the changes in total flux and heat transfer it is necessary to 
proceed to the second-order approximations, and the relevant equations for 
these are a(v? 

Vide = R pa oR sing +4 R ai oose, 
R ah ¢) | KR &R,¢)’ 
o O(9;,4,) , & (Do, 2) 
Vid, = % ye +S ae rome . (19) 
The only. difficulty in solving these equations is mechanical, and with due 


regard to the boundary conditions the second-order approximation to the 
convection is 








th, = Game 2040+ 3-75640) R2—(2-0617+-9-16070) R4+ 
+ (0°1499-+-7-5000e) R® + (1-1000 — 260000) R* — 
— (0-4250— 0-56250) R' +- (0-0342 — 0-06000) R¥® — 
— (0-0016—0-00180) R™}, 
W, = ~ iggy 08148 —8-0625 R* + 6-3406 R'—7-0583 RP + 45600 8 
— 1-7125R” + 0-3448 R™ — 0-0286 R™ + 0-0008 R'} 
+ Tama {or 154+-0-51360) R*—(0-3410+4 125210) R*+ 


+ (0-4327 + 1-14510) R’—(0-3112+.0-50000) R® + 
+(0-1312+4-0-10830) R—(0-0295 + 0-01610) R%+ 
+ (0-0024 +- 0-001 20) R'4 — (0-0000, + 0-0000,0) R'}, 
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6, = - f(0-0635 — 0-04080 + 0-13390?) — (0-1537 —0-19840+ 


~ (4608)? | 
+0:69010?) R? + (0-1914—0-39980 + 1-505907) R* 
—(0-1761 —0-43420-+- 1-82120?) R®+ (0-1103—0-27930-+4 
|. 1-34340?) R§— (0-0456 —0-11160+ 0-624207) R!4 
+ (0-0119—0-0281la+ 0-178907) R!2— (0-0018—0-0041e4 
+-0-02870?) R¥4 + (0-0001 —0-00030 + 0-002207) R' 
(0-0000, — 0-0000, 0+ 0-0000, 07) R18} 4 


, Boos 2¢ , 
T (4608)? | 





{(0-0037 + 0-05220 +-0-036607) R? — (0-0096 + 0-13180-4 


+-0:275307) R4+ (0-0106 + 0-12800 + 0-476507) R®- 
(00072 + 0-05890 + 0-369307) R8 + (0-0032 + 0-00990 4 
+- 017390) R!— (0-0009 — 0-00090 — 0-05000?) R!? 4 

+ (0-0001; — 0-00030 + 0-00820?) RM 
(0-0000, — 0-0000, o + 0-000607) R¥ + 0-0000, a? R'8}. 


These second-order approximations are included in the dashed curves of 
Fig. 3; they make no contribution to or V in vertical and horizonta! 
traverses (respectively) in Fig. 2. Flow normal to the axis of the pipe 
depends only on the parameter AB; flow parallel to the axis depends on 
AB and also on B; and the actual temperature differences across the pipe 
depend on AB, B, and c. To the corresponding order, the greatest longi- 
tudinal velocity (W) in the pipe is now at a distance below the axis 
0-6125(A B/4608)+ (0-5396+ 0-31460)(A B/4608)3 + 

Modified values of the flux of volume and of heat across any section of 
the pipe can be found now to O(A®*). If there is no secondary flow (that 
is, according to the zero-order approximation to the convection which is 


provided by the forced convection solutions) the volume flux across any 
pipe section is }vavB and the heat flux across the section (x) is 


kpc, av BiT,,(x)— rac B}, 


where c, is the specific heat at constant volume. The modified volume 
flux across any section of the pipe up to O(A?) is 


jravB{1—0-2100( 25) tees |, 
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and the modified heat flux across the section (2x) is 


rpc, av BIT, (af) —0- 2100( 55) Hap |- 


AB 
_u 5 aie 9 
brao Bl 1 (0-5177—1-02 530°)( 5) +~]}, 


from which the heat input per unit area of pipe surface is found to be 


{pe, vBr{1—0- 2100 sm) abe on |. 


The non-dimensional parameters renee of the convective flow are 
the resistance coefficient (defined as —ya/pw?,, where w,, is the mean 
longitudinal velocity across a section (cf. 1, 298)) having the value 


16, AB\* \. 
1s o4am( 4B) 


and the Nusselt numbert which is 


6/14 (0-0586—0-08520. + 0-26860?)| - a3\" fe ccohe 
| 4608 


Both the resistance coefficient and the Nusselt number are increased as 
a result of the natural convection which occurs with forced convection, 
although the increases are very small for small values of the dimensionless 
product AB of the Rayleigh and Reynolds numbers. 

Although this solution has been calculated numerically only to O[(A B)*] 
on account of the labour involved in finding higher orders, the results 
show that the forced convection solution is already in error by about 10 
per cent for AB = 3,000. Just how strong a limitation this condition will 
impose on the application of forced convection solutions to this particular 
problem may be illustrated by taking a Reynolds number 60 and a Rayleigh 
number 50. Then the following axial velocities and temperature gradients 
along the pipe wall mark a rough upper limit for the approximation (for 
this particular ratio B/A) in the case of a pipe of diameter 1 cm. For 
water, taking v = 0-01 and x = 0-0015 the greatest velocity along the axis 
is about 1-5 cm/see and the greatest temperature gradient about 0-005 

C/em. For air, taking vy = 0-15 and «x = 0-21, the corresponding values 
are about 15 cm/see and 7 °C/cm:; for mercury, taking v = 0-001 and 
x = 0-045, they are 0-15 cm/sec and 0-01 °C/em. From these very rough 
figures it is clear that the forced convection solution is inadequate even 
for quite small rates of heating of slow flows. 


+ This is defined as rate of heat transfer per unit area of wall multiplied by pipe 
diameter/k times the difference between the wall temperature and the mean fluid tempera- 
ture in a section. 
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SUMMARY 
The large-deflexion solution is given for a thin strip of lenticular parabolic section 
subjected to combined moment and torque. The analysis embraces the buckled, 
as well as unbuckled, regimes. 


1. Introduction 

Ir is well known (1) that when a thin strip is subjected to a pure moment, 
the anticlastic curvature is smaller than that predicted by small-deflexion 
theory because of the action of middle-surface stresses. Such middle- 
surface stresses arise whenever a plate deforms into a surface which is 
not developable. By the same token, when a thin strip is subjected to 
a pure torque it will assume first the (non-developable) shape predicted 
by small-deflexion theory and will later buckle (2) into a surface which 
tends to a developable one. The large-deflexion equations for a thin strip 
of constant thickness subjected to combined moment and torque have 
recently been solved by E. Reissner (3); the inherent complexity of this 
solution can be traced to the fact that the chordwise curvature of the 
distorted strip varies across the chord. 

The present paper gives the large-deflexion solution for a thin strip of 
lenticular parabolic section subjected to combined moment and torque. 
A feature of this particular cross-section is that the chordwise curvature 
of the distorted strip does not vary across the chord; this results in a 
considerable simplification in the analysis and enables the behaviour of 
the strip before and after buckling to be fully investigated. It is also shown 
that for large values of the applied moments the large-deflexion solution 
approaches that predicted by inextensional theory (4). 


2. Notation 


Ox, Oy cartesian axes, Oy measured along the length of the strip, Ox 
measured chordwise from the centre line of the strip. 
E,G Young’s modulus, shear modulus. 
v Poisson’s ratio. 
a width of strip. 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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t thickness of strip. 

D flexural rigidity per unit width. 

€ longitudinal middle surface strain, measured from a stress-free 
datum. 


w(x,y) distorted middle surface of strip. 

w(x),w chordwise variation of strip distortion defined by equation (1). 

M applied moment. 

¥ applied torque. 

K longitudinal curvature. 

twist per unit length. 

chordwise curvature (shown to be constant for the chordwise 
variation of ¢ considered). 


a strain energy per unit length due to middle-surface strains. 

V; strain energy per unit length due to flexure and torsion. 

V total strain energy per unit length. 

A constant. 

V, 2, 8,2’, M, 7,2 non-dimensional symbols defined by equation (13). 
mM defined by equation (33). 

d MT. 


Suffix 0 refers to conditions at the mid-chord (x = 0). 
Asterisk * refers to critical buckling conditions. 


3. Method of analysis 


In order to determine the relationships between bending moment, 
torque, longitudinal curvature, and twist it is first necessary to determine 
the chordwise distortion of the strip subject to a given (arbitrary) longi- 
tudinal curvature and twist. The differential equation which governs this 
distortion is determined in section 3.1. The solution of this differential 
equation for a strip of parabolic lenticular section is given in section 3.2. 
The strain energy per unit length of strip is then determined in section 3.3 
in terms of the longitudinal curvature and twist; the bending moment 
and torque may thus be readily found by differentiation of the strain 
energy with respect to the curvature and twist. 


3.1. Derivation of the differential equation 

The chordwise variation of the distortion of the strip may be determined 
most conveniently by variational methods. If the strip is subjected to 
a longitudinal curvature « and twist per unit length, @, the distortion of 
the strip is of the form 


w(x,y) = 4xy*?+-Oxry+w(z) (1) 
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where the chordwise distortion of the strip w(x), hereafter referred to 
simply as w, may be determined from energy considerations. 

The strain energy due to longitudinal strains in the middle surface of 
the strip will now be determined. From symmetry, points lying in the 
middle surface and on a given chord have the same longitudinal displace- 
ment. The longitudinal strain in the middle surface may therefore be 
obtained from equation (1), apart from a constant representing the overall 
longitudinal contraction of the strip. The first term on the right-hand 
side of equation (1) represents a developable surface and therefore con- 
tributes no middle-surface strains. The second term in equation (1) gives 
rise to longitudinal strains that vary as $6%z?. The final term in equation 
(1) gives rise to strains that vary as —«xw. The longitudinal strain 
(measured from a stress-free datum) is thus given by 


« = }$@2*—x«w+A, (2) 
where the constant A is such that there is longitudinal equilibrium, i.e. 


ta 
Bet dx = 0. (3) 


—ja 


There are no chordwise middle-surface stresses so the strain energy per 


unit length due to the middle-surface strains is given by 


ta 
Vay — 4 i} Et(}6?2? —xw+ A)? dx. (4) 


~—ja 
The strain energy per unit length due to flexure and torsion is given by 


ta 


_ (Pw(x, y) Pw(x, y) (Ae . | 
Rud | D| {Vee y)P—2(1—-»)) SF aa dx 





- ta 
ta 


+ | D{(=s) + 2vx oe ++ 2(I —vye dx (5) 


ea? 
} a 


in virtue of equation (1). 

The total strain energy per unit length is the sum of V,, and V; and the 
condition that this is a minimum with respect to w requires w to satisfy 
the differential equation 


#\p (+")|- Ext(}6%x?—xw+A) = 0. (6) 
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The four boundary conditions appropriate to this equation express the 
fact that the edges of the strip are free, whence 


d*u 
|D (G+ v| oninaiien =? | * 
d Ip d*w \ 
‘ d ae oe Yi = 0 
- last a v«)} ae 


3.2. Solution of the differential equation 


For the strip of parabolic lenticular section, we have 


f tf} -(=)" 
rnhi-C] . 


where fA en mai 5 
yp 


and it follows from equations (2), (3), and (8) that 








A —a*O?/40; 
also equation (6) becomes 


nelle] 
+ Ext} wi (= VN E— (=) | cw] —(. (9) 


[t can be readily verified that the solution of equation (9) which satisfies 
the boundary conditions (7) is given by 
Lna*{v—(Ea*6*t,/960 D,)\{}—(2a/a)*} 


w= . 5 (10) 
1+ (Ha'x*t,/960 Dy) 





and it will be noticed that the chordwise curvature «’ (= d*w/dz’) is 
constant. The longitudinal middle-surface strain is given by 
— }a?(6?+-vx?){}— (2z/a)?} 


¢€ = a re ~ . (1 1) 
1+ (Ka*x*t,/960D,) 





3.3. Strain energy in the strip 


The strain energy per unit length of strip may now be obtained by 
substituting equation (10) in equations (4) and (5). The bending moment 
M and the torque 7 may then be obtained by differentiation: 

a y Fi] 
M = ~ and JT = i= 
OK ~ 26 \. (12) 


where V=V,+), 
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At this stage, however, it is convenient to introduce the following non- 
dimensional parameters: 


0 = (FaalY =(aieek | 

~ \64 88)” ~ \(4v5)ty 

0= (“5 r) Ar a® , 
~ \(4v5)ta) 9 (assy) 


__ [(21V5)a (21v5)a 
ain (cape Ja. 64618 \r 


> a* 
“= (i) (te-ste—¢2-0 


The strain energy per unit length is now given by 
62 (82+ vR2)? 
I} 2f1+-(1—v®)R2y’ 








P= prey 


the chordwise curvature by 


ie (— —y?)§2\ 
kK = —| ————_ |k, 
1+(1—v*)r? 
and the magnitude of the middle-surface strains by 
62+-vR2 D 
1+(1—v?)a2 
All the remaining equations are deduced from these three relations. 
The bending moment and torque, from equations (12) and (13), are given 


7 m= 2 and P= (PI) (17) 


(16) 


o 
€ = 


er 2 /2 
and it will be noted that on small-deflexion theory 


M=R and T=. 


4. General relations between %, 6, M7, T 
The bending moment and torque may be expressed in terms of the 
longitudinal curvature and the twist per unit length by virtue of equations 
(14) and (17). Thus, 
7 — RI + (e+ 8%)}{1+ (1—v)(@2—8%)} (18) 
a ate {1+(1—v*)R2}2 
~ {1+ (1) 


Further, by virtue of these relations, it can be shown that 


i- ; = (5). (20) 
Ff 








and 
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It is clear from these relations and from the physical aspects of the 
problem, that for given values of 2 and @ there exist unique values for M7 
and 7, but the converse is not necessarily true (i.e. given values of M 
and 1 may be possible with more than one pair of values for 2 and §). 
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Fic. 1. Variation of curvature & with torque 7 and moment M. 
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Fic. 2. Variation of twist per unit length, 6 with torque 7 and moment M. 


This is demonstrated in Figs. 1 and 2 where curves of constant & and 
constant @ have been plotted against M and 7. 
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5. Special cases 
5.1. Pure moment 


The condition that 1 is zero and M non-zero implies that @ is zero. 
Equation (18) then yields 


uae ah ~(iaae)} enn 


The magnitude of the middle surface strains is given by 


vr 
i+ (— yr 


which tends to the value v/(1—v*) as M increases. 





“7 
€ = 


5.2. Pure torque 
The condition that M is zero and 7 non-zero implies that either 
R = 0, 
in which case f — 6{1 + (1+) 
é = & 
or 1+(1—v)(R2@—62) = 0) 
26 
l—v (24) 
l 
l—v } 


and 


in which case ? an 





and é= 


A comparison of strain energies determines which of these two states 
is the correct one. It is then found that equation (23) is applicable for 
values of |7 up to a critical value T*| where 

7* = 2(1—v)-? and 6* = (1—v)-+ (2b) 
and equation (24) is applicable for values of |? greater than \f*|. As 
the torque increases through the critical value of 7*, the strip buckles 
and there is a sudden drop in torsional rigidity. This instability is due to 
the fact that as the strip twists, the middle surface forces play an in- 
creasing part in resisting the torque and there comes a time when the 
strip will deform into a surface which approximates to a developable 
surface; for then the middle surface forces will remain constant. The 


buckled mode of deformation can be determined using the results of 
equations (15), (24), and (25): 


R2 = §2—(6*)? and # =&. (26) 


It is seen from equation (26) that & can be either positive or negative, 
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values which correspond to two distinct modes of buckling; from sym- 
metry there is no preference for either of these modes. The curvatures ® 
and ®’ increase rapidly immediately after buckling for we have 


R= +{(0—0*)(84 6%)! 
= +(S")ur *)(T 4 q*)\ (27) 








| Inextensional 
theory ~_ 
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Fic. 3. The torque-twist relationship 


Fic, 4. The torque-curvature relationship 
for pure torque. 


for pure torque. 


° . a i 7) . . . , ° 
which varies as (7'—7*)' immediately after buckling. When 7 is large 
compared with 7* A 4 

- ' Rk, Rk’ > +6 
and these relations correspond to modes of deformation which are de- 
velopable surfaces with generators at +45° to the centre line of the strip. 


The torque-twist relationship is shown in Fig. 3 and the variation of 2 
with 7 in Fig. 4. 


5.3. Bending moment and torque increasing in fixed ratio 
Consider the case in which M and 7 increase from zero in such a way 


that M Pow! 
a 
28) 
nu vam 


i.e. from equation (13) ‘= 
+y 
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The relationship between 7 and @ is found by eliminating 2 from 
equations (19) and (20): 
ge (7 /8)—1}{2—(1—v) 7/8} 
(1 +-v)[4g"{1 —(1—v) 2/8} + (2—(1—v) 2/8} 


and the variation of 2 and %’ with 7 (or M) follows immediately from 
equations (20) and (15). The relationships between 7 and 6, 2 and 2’ are 
such that for sufficiently large values of the applied moments the ratios 
7/6, 7/2, and 1/2’ tend to constant values. The asymptotic value of 7/8, 
for example, is determined by the vanishing of the factor in square brackets 
in equation (29), and is given by 


(29) 


Similarly it can be shown that 
T 2 v(1+¢?) 9 


and : m (=, Wa+9+ 08 


It will be seen from these asymptotic expressions that 
Rr’ — 82 +0 


which represents a developable surface. The angle that the generators of 
this developable surface make with the chordwise axis is given by 


ee 


which is in agreement with inextensional theory (4). It can similarly be 
shown that the stiffness, as well as the deflected shape, tends to the value 
predicted by inextensional theory as M and 7 increase. The asymptotic 
value of é is given by 


@(1—v2) > v+- {p+ y(1+4%)}-*. (31) 


5.4. Flexural snap-through buckling of a torsionally buckled strip 


Assume that the strip is in the torsionally buckled state and that the 
sign of &, determined by equation (26), is negative. The application of 
a positive moment M will then decrease |2| and at a certain critical value 
of the moment M* the strip will therefore snap through to a positive 
value of &. For purposes of calculation it is more convenient to keep 6 
constant, rather than 7, while applying the moment M. The value of 
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k* at which snap-through buckling occurs may then be obtained from the 

relation - 

OR }b const 

where M is given by equation (18) and 6 > 6*. The resulting equation 
is a cubic in &® which may be solved to give 

(1—v*)(R*)? = ph[fy(4+m)+ 24—{(4+n)—2}4]—1 ) 

where pp = {(1—v®)62—p}2 \ 

The corresponding values of 7 and M* may now be determined from 

equations (18) and (19). 


= 0, (32) 


(33) 


6. Conclusions 

A large-deflexion analysis has been presented for a thin solid strip of 
lenticular parabolic section subjected to combined moment and torque. 
The chordwise curvature of the strip is shown to be independent of position, 
and the resulting simplicity of the analysis has enabled the buckled, as 
well as the unbuckled, behaviour of the strip to be fully investigated. 
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MERIDIONAL CURVATURE 
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SUMMARY 
This paper gives a theoretical method of analysis for shells of constant thickness 
forming a surface of revolution with axially symmetrical loading and with constant 
meridional curvature. The method of solution is to solve directly the basic 
differential equation governing the problem. The method is comprehensive and 
applies to a wide range of design problems, although in the present instance it is 
used to determine the stresses in pressurized corrugated ducting. 


1. List of symbols 


M, = Meridional bending moment per unit length. 

M, = Circumferential bending moment per unit length. 

Ns = Meridional force in shell per unit length. 

No = Circumferential force in shell per unit length. 
€g = Circumferential strain. 

Q,4 = Shear force per unit length. 
r, = Radius of curvature in meridional plane. 

r. = Radius of curvature in the normal plane perpendicular to the 

meridian. 

r, — Radius of shell. 
R = Maximum radius of shell. 
¢ = Angle between normal to shell and shell axis, 
pb = $n—4¢. 
p = Pressure. 
h = Thickness of shell. 

E = Modulus of elasticity. 
v = Poisson’s ratio. 

D = Flexural rigidity of shell. 
V = Angle of rotation of a tangent to the meridian. 
U= Q% T>. 

F(p) = Pressure function (equation 10). 
L{...) = Differential operator (equation 9). 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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p = [ £h/D—v?/r?}t. 
S,, S,, S3, 8, Power series in ~. 
A,, Ay, Az, A, Constants of integration. 


n = Solution of equation (26). 
y = re(n). 
5 = im(n). 

71> Nz See equations (34) to (36). 
r = Positive integer. 

a, = Coefficient in power series. 


° 


= Constant in power series. 
K = rir, p/2nD. 

A = —(5+2v)r3 p/2nD. 

n = 7? u?/ro. 


2. Introduction 

The mathematical theory for cylindrical shells and spherical shells is 
comprehensive and, with a few specific exceptions, adequate. However, 
the theory for shells of other shapes is by no means as complete. In 
particular, shells forming a surface of revolution with constant meridional 
curvature, such as toroidal shells and corrugated shells, are not covered 
very adequately by existing theories, although recently work has been 
done by Clark (1) and by Turner and Ford (2). The former work consists 
of a comprehensive analysis which uses asymptotic integration, and conse- 
quently the method is not valid for the same range of parameters as the 
present method; i.e. the present method would apply to small corrugations 
whereas the method of Clark would apply to relatively large corruga- 
tions. The latter is an approximate technique using minimum strain 
energy, whereas the present method uses a direct approach to solve the 
relevant differential equation governing the problem. The method of 
Turner and Ford is applied to the corrugations occurring in bellows units 
where the semi-included angle of the corrugation is quite large, and may 
be greater than 47, whereas the present method is developed for corruga- 
tions occurring in flexible pipes and ducting where the semi-included angle 
is not likely to exceed }m or 47. 

With such problems, where bending in the shell is important, a solution 
whereby the membrane stresses are first obtained and then the bending 
stresses are added is not justifiable, since the deflexions implicit in the 
membrane stresses imply appreciable bending. This point has been clearly 
brought out by Fessler and Rose (3) for an elliptical shell, although it has 
been appreciated by many writers. 
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The present theory is concerned with axial symmetrical loading and 
the starting-point of the solution is the equilibrium equations for a shell 
element, that is in effect, the Love—Meissner equations. The solution 
does not require the separate determination of a membrane solution, and 
the resulting differential equations are solved in terms of series. The 
mathematics of the problem is somewhat similar to the case of a corru- 
gated diaphragm which has been solved in great detail by Stange (4). 

The work is done principally in order to investigate corrugated ducting, 
since little work has been carried out on this particular subject. 


3. The general mathematical problem 

For any shell forming a surface of revolution and loaded by means of 
internal pressure, the equilibrium of any small element leads to the 
following three equations (see Fig. 1 and (5)). 


CIRCUMFERENTIAL SECTION. MERIDIONAL SECTION. 











AXIS OF SYMMETRY. 





Fic. 1. Element of shell. 


dw 
dg Noro) —Nors 008 $—10 Oy = 0, (1) 
Norot Nori sing +7 (Qpre)—rarop = 0, (2) 
d 
dg Mora) — Mor, cos 6—Qg 1170 = 0. (3) 


The longitudinal force in the shell, Nz, can be expressed in terms of 
Q, and p, either by elimination between equations (1) and (2) or, since 
Ng is constant along a parallel circle defined by the angle ¢, by a considera- 
tion of the equilibrium of that portion of the shell lying on one side of 
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the parallel circle. Hence, 


= —Q, cot ¢+—”" (4) 


2 ied 
It should be noted that for the shell under consideration the pressure 


load on the shell ends is reacted by the shell itself. That is, the total 
axial load on the shell is prr?. A sketch of the shell is given in Fig. 2. 


"i. 


ee 








AXIS OF SYMMETRY 





Ta Mal ia eas 
6 


| 
| 
} 


____axts'OF SYMMETRY. _ 
Fic. 2. Type of shell. 


Other loading cases can be analysed in a similar manner, the minor 
alterations required being outlined in section 3.3. 
If the expression (4) is substituted in (2), Ng can then be found: 


By = > 5 (Qyry)— Oh 4 0 (5) 


2r,sin’d ' sind 


+ 


If the shell is of constant thickness, then, by (5), the rotation of a 
tangent to the meridian is given by 


Ehr,V = cot d[(r,+r.)Ng—(r2+ vr )Nol— ae! ‘o(Ng—vNg)] (6) 
and by substituting (4) and (5) into (6) we find 

rs d*(Q4 Ts) +0 d(Q472) Lir, a 

Bae te dagle) tote |g 7 [peo] or 


“ env +| 7 2 (5-+2v)r Jae (7) 
1 





or L(U)+=U = EhV+F(p), (8) 
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we: rgd...) Lf d (re\ d(...) cot®¢(...) - 
9 tS tate ew 


ip) = [78 (6+ 29g] PONE: (10) 


the second term in F(p) is smaller than the first, but not negligible in the 
present instance. Similarly, if 


r, dp 


V outg v | 
‘2 Th d¢ 


are substituted into the third equilibrium equation, then the resulting 
differential equation is 


M, = —D|- +2 Veot¢| (11) 


and My = ak (12) 


, U 
1V)—~V = -> (13) 
The two simultaneous differential equations (8) and (13) are as derived 
by Timoshenko (5) and others, for the case p= 0. However, before 
considering the method of solution of these equations, the differential 
operator L(...) can be simplified for the present problem. In particular 
for the case under consideration, r, is constant but r, varies; in fact, 


r,8ing = fp 
and ro = R—r,+r,sin¢ 
where R is constant. Hence 


ti d*(...)  cosd d(...) _cos*¢(...) 
HA) om r?sing dd? TF sing df = rasing ” (15) 
where it should be noted that the coefficients are of a reducing magnitude 
provided r,>r,. Although not necessary for the validity of the solution 
the numerical work can be simplified by using the approximation 


rt, d*...) 
rising dg? 
which is justifiable provided rz > 7,. Further comments are made abo: t 
it later. | 

The second-order simultaneous differential equations (8) and (13) can 





Li...) = (16) 
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be used to give fourth-order differential equations in U or V. These 
are 


ae i. ‘ ; 
L.Lt + | 5A! = UF(p)) —" FP) (17) 

ase ee 
and L.UV)+| 5 —4|¥ = —5 Flo). (18) 


To solve (18) we require the complementary function, a solution of 
L({V)+p4V = 0 (19) 
together with a particular integral of 
LV) +pV = ~7 F(p). (20) 
3.1. The complementary function 


The solution of equation (19) is well established; the equation breaks 
up into the equations 


L(V)+142V = “ (21) 
and L(V)—ip2V = 0 r 
where if V =8,+18, | (22) 
and V =8,+%S, j 
(S,, S,, Ss, and S, being real functions) are solutions of the first of equations 

» y Y 2 

(21), then V =S,—i8,) (2) 
and V = S8,— iS, 


are solutions of the second of equations (21), and consequently the com- 
plete solution of equation (19) is 


V = A,S,+A,8,+A,;8,+4,5,, (24) 
where A,, Az, As, and A, are real arbitrary constants. 


3.2. The particular integral 


The particular integral of equation (20) can be derived by reference to 
the work of Stange (4). Consider the complex function 


n= y+, (25) 
which is a solution of 
‘ l 
n= — ; 26 
L(y) +n 2D F(p) (26) 
By considering real and imaginary parts, we find 
l 
Ly)—p = — 
(y)— 2D F(p) (27) 


and L(8)+p*y = 0 
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and hence D*(5)+p% = — 5 Fv). (28) 


Hence V = 4 is a particular integral of equation (20), 5 itself being found 
by solving equation (26) for ». 
3.3. Other load conditions 


For other load conditions only minor modifications are necessary. For 
example, a pure end load of magnitude L, would give 


Ly 


2rr,sin d 


Le 
2er Dar, sin’ 


from which F(p) = E —(3+2y on 
Similarly for an internal pressure with no end load, 
cot ¢' 
F(p) = — (roe 


7 


= —Q,cotd+- 


and No = <2 ~ 14 (Q- -T.)— 
1 


2r3 


3.4. Displacements 
Displacements can be obtained, if required, by integration of the ex- 
pressions for strain, i. 
T1€4 = —-— W, 


dd 
that is, r€g = veotd—w. 
The resultant longitudinal displacement v, and radial displacement w, 
can be simply deduced from the solution of the above equations. 


4. The method of solution 

The basic differential equation is linear with variable coefficients, 
and the solution can be obtained by means of infinite series. The main 
advantage of this technique is that it is general and is likely to be 
practical for a wide range of functional forms for both the differential 
equation coefficients and the particular integral terms, F(p). 

The paper by Stange (4) provides a useful reference since, whilst the 
physical problem with which he is concerned is somewhat different, the 
mathematics of both problems is similar. 

For the complementary function the equation to be solved is the first 
of equations (21), and the approximation given by equation (16) is used. 
As mentioned earlier, the method of solution used here will still be 
applicable if this approximation is not made. The advantage of in- 
cluding it is that for very little loss in accuracy there is a large saving in 
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the numerical work. However, it should be remembered that if r,/r, is not 
negligible compared with unity, then the full differential equation should 
be used and the technique employed here will still be applicable. 

For the particular integral, a similar technique is used in order to solve 
(26). 

4.1. The series for the complementary function 

Substituting (16) in the first of equations (21), the differential equation 
for solution becomes 

=e . ae 
2ainddd2 
ri sin d dd 

{ys 2, 23 

4 ae ain¢ sy — 9, 
dd ro 
and by putting n = rj?/r, (where r, can now be given its mean value in 
the case of the corrugated ducting under consideration since the variation 
is small) and ¢ = $xr—4¢, 


27) = 0, 


i.e. 


@2V 
di? 


Further, as a close approximation, 


+inVecosy = 0. (29) 


yy 


cosy = 1— ait 4)” 


this in fact being true to 1-4 per cent. for 4 = 1. Thus 


2V 2 4 
pe +in(1 5th)? =’ (30) 


Now V can be expanded as a power series in ys, which is certainly convergent 
for the present range of yf, i.e. 0 to approximately 1. Thus, putting 


V few: 2" (31) 


and substituting in (30) it is found that c = 0 and that all the coefficients 

of even powers of % can be expressed in terms of ay, and that all the coeffi- 

cients of odd powers of y% can be expressed in terms of a,. This leads to 

two complex series in powers of y. If ag = a, = 1, these series are 
S,+i8, 

and S,+7i8, 

so that as explained earlier, the solution for V in equation (19) becomes 

V = A,S,+A,S,+4,5;+4,S,, 
where the right-hand side is a real series in positive powers of #. 
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The actual valics of the first few terms in the series are as given below 


8, = 1- ni 4 int + (nt—3int) © + (—50n'+ 98n2) 


10! 
_ +, ¥ so 3¥* 73) 0° 
Ss, = —n aj + gy t (mm) Gi — 2208 + (0+ 2078!) 
_¥ 


3 
S 


5 7 9 
no + lant + (n*— 1038) & + (—70n*+ 558n?) v 
Vv. 3 : 


11 
Oe 
+ 3a 


uu 
! 





yp yf _ yi 
+ (n3 . la 37 4 (m5 § 
(n?—5n) 71 34n oi! n>+-697n°) im | 
4.2. The series for the particular integral 


For the particular integral we require any solution of (26), this equation 
being rewritten as 





dy. _ og Sing sin ys 
aga t m7 008 = 2K ee ot 
Putting 7 = (+ tan p)K + (n_+log{tan(}6+ }7)})A, (34) 
the resulting equations for yn, and yn, become, approximately, 

dm. a a PO 

ai + inn (1—5, +4) ca ~in(y 5 +5) (35) 
d ; 
and Tg tinny) 7 =}: (36) 
These two equations can in turn be solved by series and the constants 


can be chosen arbitrarily, since only one particular integral of each is 
required. The first few terms in the series for , and , then become 


(33) 


5 7 9 
= | nti + lint + (—s7mt tn ST 4 


a ee ee y 
+i[ mF +m +( ntnty ts +(n—32n)5 | 


3! 





ms = —nt 4 ant | +i an +20 +m] 


5! 
5. The complete solution 
The complete solution for V is 

V = A,S,+A,8,+A,8,+4,8,+5 (38) 
where A, to A, are arbitrary constants, S, to S, are the series given by 
equation (32), and 8 is given by 

7 = y+ (39) 

where 7 is given by (34) and (37). 
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The complete solution for U need not be evaluated separately, but can 
be deduced from that for V, by using equations (13) and (21). In fact, 
if this is done, then the complete solution for U becomes 

! Y + Agu®)S,+(4,” —A,p*|S 

a A,—+Agu wd Eel Set 

1 1 
iy (4, ~ + Aun’) +(4e ~ —Agp8) Stuy + ~8. (40) 
" r " 

Thus if V and U are completely known, then M,. Mg, Ny, No, «4, and «¢ 
can be evaluated, since the latter six functions can be expressed in terms 
of the former two functions. 


6. The boundary conditions 

The complete solution of the differential equation of the problem 
involves four unknown constants, so that four boundary conditions are 
necessary to specify the solution completely. 

For the case of a corrugated pipe, it can be seen, in Fig. 2, that the 
solution is applied to two sectors numbered I and II. In sector I, 7, is 
positive and in sector IT, it is negative. 

For the case under consideration the boundary conditions are as follows. 
Firstly, by symmetry, at ¢ = 0, V = U = 0. With reference to Fig. 2, 
this applies to both sections I and II. Secondly, by continuity at 4 = y,,,, 
V and U for section I must be equal to V and U, respectively, for section 
Il. Also My, and «9, or if more convenient, two other functions can be 
equated at this junction. 

These conditions lead to the evaluation of the A’s, and in fact, by the 
symmetry, A, = A, = 0, which results in a simplification of the numerical 
work, 

For other problems for which the solution is applicable the boundary 
conditions can be applied in a similar manner. 


7. The numerical results 


The complete solution has been evaluated numerically for the case of 
the particular corrugated duct for which 


ro = 17-5 in. (mean)) 

r, = 1-791 in. 

100 Ib/sq. in. }- 
h = 0-375 in. 

vs, = 0-993 rad. } 


The bending stresses and the direct stresses have been found and com- 


+ 
! 
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| CIRCUMFERENTIAL STRESS. 
Fics. 3 (a) and (6). 




















bined; the resultant longitudinal and hoop stresses on the inner and outer 
surfaces are plotted in Figs. 3(a) and 3 (6). 

The above numerical values are those which correspond to an actual 
corrugation of 36 in. nominal diameter, axial distance peak to trough 
3 in., radial distance peak to trough 1§ in., and inside radius of meridional 
curve 1} in. 

5092.48 Gg 
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This particular corrugation has been tested and the stresses measured 
experimentally, the experimental values being also plotted on Fig. 3. 
The circles giving the stresses on the inner surface and the squares being 
the stresses on the outer surface. 

This shows quite satisfactory agreement, particularly when it is remem- 
bered that the calculations apply to mean properties of the corrugation, 
whereas these can be expected in practice to vary throughout the speci- 
men—the stresses, in particular, being very sensitive to wall thickness. 
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SUMMARY 

The behaviour, under transverse load, of a flat plate stiffened by beams is analysed 
by Fourier series and solutions are obtained for symmetrical and antisymmetrical 
forms of loading. 

An alternative method of analysis is derived for symmetrical loading forms, in- 
volving the evaluation of the plate efficiency. This alternative method has some 
advantage in that the applied loading is expressed in the form of a Fourier series 
for bending moment rather than load intensity. 

Finally, it is shown that approximate solutions (assuming infinite torsional rigidity 
of the beams) may be obtained for a wide variety of loading forms by superimposing 
the solutions for the symmetrical and antisymmetrical forms of loading. 


1. Introduction 


THE analysis of ‘deck systems’, either in the form of interconnected beams 
or stiffened plates, is a complex problem which has been under investiga- 
tion for some considerable time. Recent advances in the theories relating 


to the elastic behaviour of interconnected beams have been described (1, 2). 
This paper considers the related problem of the stiffened plate, involving 
the shear lag effect and plate efficiency. 

When a stiffened plate structure is subjected to transverse loads, shear 
deformation occurs in the plate which results in a variation in longitudinal 
stress across the width of the plate. The longitudinal stress in the plate 
is a Maximum at its junction with a stiffening beam. This phenomenon 
is generally termed the ‘shear lag’ effect. 

Under symmetrical forms of applied loading the maximum longitudinal 
stresses and deflexions occurring in a stiffened plate structure may be 
determined by considering the ‘effective width’ of the plate (3, 4). Fig. 
1 (a) illustrates a semi-infinite, stiffened plate under symmetrical trans- 
verse load; thus the stiffening beams are identical and equally spaced, and 
are each subjected to the same applied bending moment. For this form 
of structure and loading a ‘repeating element’ may be isolated which 
consists of a stiffening beam and a width of plating equal to the beam 
spacing. The efficiency 7 of the plating in forming a flange to the stiffening 
beam may be defined as the ratio of the average to the maximum longi- 
tudinal stress in the plate; the effective width of the plate is then given 
by the product of the efficiency » and the stiffening beam spacing. 


[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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When analysed by the elementary theory of bending, the maximum 
longitudinal stresses and deflexions of the composite beam (formed by a 
stiffening beam and an effective width of plate) are the same as those 
which occur in the actual stiffened plate structure. For forms of loading 
which are not symmetrical and do not give rise to a repeating element 
(e.g. see Fig. 1 (5), (c), and (d)) the plate efficiency analysis does not apply. 
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In this paper an expression is derived for the plate efficiency » applicable 
to the type of loading shown in Fig. 1 (a). The analysis given does not 
depend on the assumption of widely spaced stiffening beams or relatively 
thin plating, as do previous methods which have been developed (5, 6). 

Complete solutions, as distinct from the more limited solution for plate 
efficiency, are subsequently derived for the symmetrical and anti-sym- 
metrical loading forms illustrated in Fig. 1 (a) and (b). Finally, it is shown 
that, by superposition of these type solutions (7), approximate solutions 
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may be obtained for many other loading forms (e.g. those shown in Fig. 
1 (c) and (d)). 


2. Notation 


The form of load applied transversely to the stiffening 
beams p = > p, sin(n7z/l). 
The form of bending moment applied to the stiffening beam 
M = > M, sin(n7z/l). 
Modulus of elasticity. 
Poisson’s ratio. 
; Modulus of shear. 

Sabates Plane stresses in the plate. 

u,v, w Displacements in the z, y, z directions respectively. 

M,, M,, M,, Bending moments per unit width in the plate. 
Reactions per unit width in the plate. 
Moment, normal force, and shear force carried by the stiffen- 
ing beam. 

F,(F,,,) Average longitudinal stress in the top surface of the plate 
due to the given loading (due to harmonic component). 

F,,( Fn.) Maximum longitudinal stress in the top surface of the plate 
due to the given loading (due to harmonic component). 
Effective width under the given loading. 
Effective width for a single harmonic component of loading. 
Depth to neutral axis associated with flange width 7, a. 
Second moment of area associated with flange width 7, a. 
Span of stiffened plating. 
Spacing of stiffening beams. 
Thickness of the plate. 
Distance from centroid of the stiffening beam to centre of 
plate. 
Cross-sectional area and second moment of area of stiffening 
beam alone. 
Dimensions of a rectangular cross-section stiffening beam. 

D Flexural rigidity per unit width of the plate. 


3. Stress distribution and displacements in stiffened plating 


The general notation and sign convention for the analysis of the sym- 
metrical loading case of Fig. 1 (a) is shown in Fig. 2. An elemental length 
of the stiffening beam is shown in the figure, indicating the forces 
and moments exerted by the plate on the beam element. The suffix ‘0’ 
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indicates that the quantity concerned is taken along the line of the beam, 


i.e. at y = +4a. 
The stress components f,,f,, and f,,, must satisfy the equation 
VF = 0, (1) 
A solution for f, which satisfies equation (1) is 
~ nry pry ny nry 
A, cosh - . A, sinh +A, cosh —= + 
Se doa 1 i sj 2 l l 
n= 1,3,5.,... 
+A —— sinh 74 \sin"™~ — (2) 
l l 
in which A,,..., A, are constants to be determined. The plate stresses /,, /,,, 
and f,,, must also satisfy the following conditions of strain compatibility : 
of, A fry 0 ) 
Ox ly 
Sev 5 Fv _ 9 
Ox oy 
ou 
E = f.—; , (3) 
‘tad ay ie | 
as 
E= re 
cy | A f, 
Ou ov 
Patel Ha 
lay) | 





Symmetry about the axisy = 0 gives A, = A, = 0, and using equations 


(2) and (3) and the conditions v = 0 at y = +4a and 0 gives 
A, Zi 3+v nna @ oth U7. (4) 
A, i+» WD 21 


The displacement of the plate perpendicular to its own plane must 
satisfy Thy = @, (5) 
where second-order effects due to membrane actions have been neglected, 
and for simplicity the applied loading is considered concentrated at the 
stiffening beams. A solution for w which satisfies this equation is: 


x 


© = (2: cosh 4 4-B, sinh e+ Bs re osh 74 4 
n=1,3,5,... l 
B, ney nh MTY) ain ™™™ 16 
+B," Y) sin™™™ (6) 


in which B,,..., B, are constants to be determined from the boundary 


conditions and nm has odd values only. As w must be symmetrical with 
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respect to the axis y= 0, it follows that B, = B,=0. The further 
condition that éw/éy = 0 at y= +4a and 0 gives 
z = -(1 +" coth (7) 
An additional relationship may be established by equating the strain 
in the plate and stiffening beam at their common surface. From inspec- 
tion of Fig. 2 this relationship may be expressed by the equation 


a2 . , 
Hw, , N 


: Cw, 1 ; 
(9) 0 + Fy = 0 Ef, 


Ser Ye ax 






































Fic. 2 
The individual terms in equation (8) may be evaluated from equations 
(2) and (6) and the values of the coastants A,,..., A, and B,,..., B, so far 
determined; in addition by considering equilibrium in the z-direction we 
find | eee [Savy de. 
The solution of equation (8) leads to the result 


BR, tte" 
7 Ads) Be 
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hens Leis S—o— (1-0 A+ {ahi Ann(1 +7) eonh(nwe/2) 
(1+A)e/d 








nra/2l 


- 7 sinh(n7za/2l)cosh(n7a/2l) 





Thus the eight constants appearing in equations (2) and (6) may all be 
expressed in terms of one unknown, say A,. 
The plate efficiency may be expressed by the equation 


F. 
| F,, 
where als 
l 6M 
m7 *= [ (— = +h, dy 
ee 
“ ., nnalh{ A 21 4l . nx 
= A.t vereciace! fe n_\{_: aa Pne/ae 
>, os 21 lin tallona) tac a mK l . (10) 
and i. = -..™ +f, 
” h2 
= n h\ A 
= bY A, cosh 7 (aaa ty 4 A)+2v}+ 
nodd - 
3-+y . nrx 
+ —— —A|sin —_— 
l+yv é 4 





From equations (10) the plate efficiency under a single harmonic com- 
ponent of bending moment series may be shown to be 


h 4 l na 
> hee a os - |tanh — 
lalsci-2-n) mel m3 


h A , 3-+y 
——™ —_ {(] —y)(1+- +2v}4 ai 
d xi — a)" r(L-+A)+ 2} l+y a| 





Mn (11) 


na 


21 





as the unknown constant A, disappears. The plate efficiency under an 
applied bending moment of the general form M = > M, sin(nzz/l) may 
thus be determined by superimposing the effect of each harmonic com- 
ponent of the bending moment series. The process is briefly described 


below. Thus - F,, zs > F,, % > F, (12) 
U] — F — > F — ye : ° - 


To ne 


For a known distribution of applied bending moment the functional form 
of M,, may be determined and sufficient terms taken (e.g. n = 1, 3, 5,...) 
to give an accurate series representation of the applied bending moment. 
Values of », may then be determined for each value of n from equation 
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(11); for each value of »,, the corresponding values of d,, and J, are found, 
and thus, using the elementary theory of bending, 


F,, { = (d,, M,,/I,)sin(nz2/1)] 


may be calculated for each value of n. Finally, the plate efficiency may 
he found from equation (12). For a given applied bending moment the 


Sr 


es 
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plate efficiency may be calculated in the above manner for any value of 
x from 0 to /. In a similar manner the stiffening beam defiexions may be 
determined thus: 


l l \?M, si l 
~o=Tu=5 (--) aaa) 


In Fig. 3 values of » are given for two forms of applied bending moment 
involving simply supported spans. If restraint against rotation of the 
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structure exists at the supports, solutions may be obtained by super- 


imposing on to the simply supported solutions the effects of a constant 
bending moment. 


4. Complete solution for the symmetrical loading form of Fig. 1 (a) 

The analysis for determining the plate efficiency may be extended to 
give a complete solution for equations (1) and (5) by evaluating the 
constant A,. This constant may be determined by considering the equi- 
librium of the beam shown in Fig. 2, and three equilibrium equations may 
be obtained thus: 


In the z-direction pdx+dT—2T, dx = 0. 
In the x-direction dN +-2hf,,,dx = 0. (13) 
Moments about the beam centroid 
dM’ = Tdx+-2hef,,,dx = 0, 
d*wy 
dx? 
Elimination of M’, NV, and 7' from equations (13) gives the result 


where M' = —E!I 


OMe os & 
J) = 2h ~ rVo — > _ = 0 
El dat he Te in p ' (14) 


where p is the applied load intensity expressed in the general form 

p = > p, sin(nzz/l). 
The individual terms in equation (14) may be evaluated from equations 
(2) and (6) and the values of the constants A,,..., A, and B,...., B, so far 


determined. The solution of equation (14) then gives the following value 
for A,: 


IB 
A ee la Pn ee : 15 
. nrd* cosh(nza/2l) —" 
where 
I (nza , 4(h/d)* tanh(nza/2/) 4he tanh(nza/2l)]-! 
apes Pe 14a)4 34 ce as med sate: Lew | 
“ alaral l | 1—y? d*(1+-v) | 


The basic equations (1) and (5) which define the stress distribution and 
deflected form of the stiffened plate are thus fully solved. In applying 
this solution it must be ensured that convergent series are used to express 
applied load intensities in the form p = > p, sin(nzz/l). In particular, 
series representations of point loads and rotational restraints at the sup- 
ports give difficulties which are not without physical significance. A point 
load, in this connexion, cannot truly exist but must be distributed over 
a small finite distance. The Fourier series for the load intensity shown in 
Fig. 3 represent high-order load intensity distributed over a very small 
fraction of the span, rather than point loads where the load intensity would 
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tend to infinity under the load. The plate efficiency analysis does not 
suffer from these difficulties as fundamentally the process is one of super- 
posing the effect of each harmonic component of a bending moment series. 


5. Solution of the antisymmetrical loading form of Fig. 1 (6) 
The method of solution for this loading form is similar to that used for 
the symmetrical loading and will, therefore, be only briefly described. 
Solutions for equations (1) and (5) are, respectively, 
A \ 
f. > (c; cosh ” “¥ +0, sinh —t+0, met cosh “TY + 


nodd 


+0,-2 sinh "7s in “= 


. (16) 


> D, cosh“ + D,sinh “4D, mrt cosh TY + 
nodd 
+" sinh “sin 
Due to antisymmetry about the axis y = 0; C, = Cy = D, = D,= 9; 
v = Oat y = +4a; and éw/éy = 0 at y = +4a. Hence 





Cy we OP nt * tanh ("se 
ry l+v 2 2l 


(17) 
and 7 -(14+55 re. ~. oT “ta anh") 
The relationships established in equations (8) and (14) still apply, and 
yield the following values for D,/C, ae Cs: 
Dy _ 2D, \ 
i Ed 
ists peat + es piais ie Ae asee 
“ (1—A)e/d 
E, lp, . (18) 


~~ eosh(nza/2l) nad? 





and 


where 
. {i na na (2 -1 
E, - | ile (1—A)tanh | _ pay | 
6. Superposition of the symmetrical and antisymmetrical loading 
forms 
Figs. l(a) and 1(b) represent the two loading forms already fully 
analysed. Let the effects (including deflexion, stress, etc.) produced at 
the stiffening beams for these forms of loading be s and ++ respectively, 
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where s and ¢ may be calculated from the foregoing analysis. Superposition 
of these loading forms gives the loading condition shown in Fig. 1 (d), 
where the effects will accordingly be }(s+-t) and }(s—t) at the loaded and 
unloaded beams respectively. 
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Now consider Fig. 1 (c) where only one beam is loaded, and let the effects 
occurring at successive beams outward from the loaded one be «, f, y, and 
zero. By reconsidering the effects under the loading shown in Fig. 1 (d) 
in terms of a, 8, and y, and assuming infinite torsional rigidity of the 
beams, the following equations are obtained: 


28 = 4(s—t), a+2y = (s+), (19) 
therefore B= }(s—b), x = $(s+t)—2y. 
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Hence 8 may be determined directly in terms of s and t, In the majority 
of real structures a> y and, therefore, an approximate, conservative 
value of « is $(s+t). In this way the symmetrical and antisymmetrical 
solutions may be superposed to give approximate solutions for a variety 
of loading forms. 


7. Experimental results 

The results of tests carried out on several small-scale stiffened plate 
structures confirmed the validity of the above analysis and the assump- 
tions on which it is based. Fig. 4 shows some of the experimental results 
from one of the specimens. This test specimen was constructed in 
aluminium alloy, with the following dimensions: a/l = }, h/d = 4, b/d = 4. 
In Fig. 4 the experimental points and corresponding theoretical relation- 
ships (full lines) are shown for concentrated loads W applied at the mid- 
span of the simply supported specimen. 
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SUMMARY 

It is shown in this paper that simple assumptions about the transverse displace- 
ment and the shear stresses lead to two equations for flexural vibrations of plates, 
corresponding to three modes of vibration. The equations obtained are compared 
with the exact equations for the propagation of plane flexural waves, and it is 
shown that if the frequency is less than a certain value, correspondence is very good 
between the two sets of equations. The results are also compared with similar 
results obtained from Mindlin’s equations. 


1. Introduction 

THE aim of this paper is to obtain an approximate equation for flexural 
vibrations of an elastic plate for fairly low frequencies and to find out 
just what assumptions have to be made about the stresses and displace- 
ments in order to obtain a reasonable theory. The fundamental approxi- 
mation of this paper is that the shear stresses vary parabolically across 
the thickness of the plate, an assumption which is certainly correct for 
very long waves, as has been demonstrated by Prescott (1). Reissner (2) 
makes a similar assumption and the theory developed in this paper is, 
therefore, analogous to Reissner’s theory of static bending, although the 
method used to derive the equations is rather similar to A. E. Green’s 
derivation of Reissner’s theory (3). The notation, however, is that of 
Stevenson (4), and it has been found necessary to make a further assump- 
tion about the transverse displacement, which is expressed in the form 
of the first two terms of a Taylor series. 

The accuracy of any theory of plate vibrations can be checked by 
comparing it with the exact equations for the propagation of waves in 
plates. The equations obtained in this paper are found to agree very well 
with the exact theory for fairly low frequencies, although not significantly 
more so than Mindlin’s theory of plate vibration (5). Mindlin’s theory, 
however, is based on unrealistic assumptions about the displacements, 
assumptions which are compensated for by an undetermined constant, the 
value of which is fixed only after comparison with the exact equations. 
Whilst it is not suggested that the equations obtained in this paper 
should supplant Mindlin’s equations, they are certainly based on a better 
foundation. 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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2. The fundamental equations 

It is assumed that the mid-plane of the plate is the plane Z = 0 of 
a set of rectangular cartesian coordinates (x,y, Z); we denote by u, v, w 
the displacements in the directions O(z,y, Z), and the faces of the plate, 
of thickness 2h, are the planes Z = A. 

The complex stress combinations ¥’, ®, © are defined by 


Y = 2z+iy2, (2.1) 
O = zr+ 99, (2.2) 
and © = rx—yy+-2izy, (2.3) 
whilst the displacen-ent D is defined as 
D = u+w. (2.4) 
The complex stress resultant % and stress couples I’, Q are defined by 


(2.5) 
(2.6) 


h 
and (2.7) 


(2.8) 


ow 


éD 
= p— —» 2. 
¥ ry ia (2.9) 


(1—20)@—2022 = 2u(2 +2) 


(2.10) 


and (1—o)z2—ocO — (2.11) 
whilst the equations of motion are 


oo © oF ak 


tate = (2.12) 


oY oF az s 
and me tot eZ a pw, (2.13) 


where p is the density, o = A/(2u+3A), A, w are Lamé’s constants, 
z = x+iy and the dots indicate derivatives with respect to the time t. 
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The boundary conditions on the faces Z = +h are 
2a = 2y = 22 = 0. 
3. The derivation of the approximate equation 
It is assumed that Y = 34(1—Z?/h?). (3.1) 


The assumption that w is constant with respect to Z is insufficient for 


this theory, and Reissner’s average w* leads to insoluble equations. We 
assume that 


Ww = Wet }w,(1—Z?/h?), (3.2) 
where w, and w, are functions of 2, y, and t. Integrating (2.12) with 
respect to Z between the limits =h, and using (3.2), we get 


+e. 


= = p(t) +t). (3.3) 


Substituting the expressions for ‘Y and w from (3.1) and (3.2) into (2.13), 
we find, with the use of (3.3) and the condition zz = 0 on Z = $A, that 


22 = tpi, Z(Z*/h?@—1). (3.4) 
Multiplying (2.9) by (1— Z?/h?) and integrating we obtain 


uy 5uD* | 5 
¥= “oe + 3M 


po + ay =. 
h ¥ 
where D* = - [ DZ dZ. (3.6) 


If (2.8), (2.10), (2.11), and (2.12) are multiplied by Z and integrated, it 
is found that aD* 


lr = 4n-—_, (3.7) 
0 
he bi lle = 2u(1+n)é, (3.8) 
and a. (3.10) 
Cz 

where » is Poisson’s ratio and 
(3.11) 

From (3.8), (3.10), and (3.11), 


a o 
uv2p* +2," +7) of 2pnh as = pb*, (3.12) 
l—n @ 15(1—n) Oz 


—_ a 
Toya 


where V? is the two-dimensional Laplacian operator — Separating 
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the imaginary and real parts of pr we obtain the two equations 


co * tay 
where (= wot _ out D* = u*+iv*, 
ox hy 


€ 2 es 

and = sl Vig —p(tig t+ ,) = p&. (3.14) 
l—y 15(1—7) 

Taking the real part of (3.5) and using (3.3), 


pltig+té,) = seat + §uV2w + pV wy. (3.15) 


Equations (3.9), (3.14), and (3.15) form a system of three equations for 
Wo, W,, and &. If these three quantities are expressed in terms of harmonic 
components, it can be shown that if w, = O(1), then w, and é are O(h?/A*), 
where A is the wavelength, and that the operators h?V? and @*/ét* are also 
O(h?/A®). Retaining only terms up to O(A®/A®), (3.9), (3.14), and (3.15) can 
be reduced to the single equation 
2uh? ; 17—7y _ 67? +58n—59 lv) “~ie 
3(1 om bbe ~ pa 7? ph?V*i, 150(1—y)e 2*h*wi!*) + pti, = 0. 
(3.16) 
If c, = (u/p)' is the velocity of shear waves in the material making up 
the plate, and c, the velocity of compression waves, with y = c?/c3, then 
(3.16) can be written in the form 
4(y— D 2V4— oTy— 20 sevay + 57y?—50y—10 giv 4 Wo = @, 
By L5yc3 150y(y— 1)ec$ 2 





(3.17) 
whilst (3.13) may be written as 


5 ¢ 
on ee = a (3.18) 
Equations (3.17) and (3.18) are the equations for flexural vibrations, [ 
being independent of w, unless coupled by the plate boundary conditions. 
¢ consists purely of rotations about axes parallel to the axis of Z, whilst 
w, and € are displacements consisting of the dilatation and rotations 
about axes parallel to the plane of the plate. 


4. Comparison with exact theories 
If we put Wy = wy etvz-m, (4.1) 
where wy, is a constant, then we obtain the following approximate wave 
equation for the propagation of fiexural waves in elastic plates 
2 hale 1 — 
wt 4(y—V) pape 27Y—20) pap ays + (iv — kent, (4.2) 


3y L5yc} 150y(y— 
Hh 
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where the wavelength A = 27/p, and n/2z is the frequency, whilst if 

t = C’etpr- nb) (4.3) 
the wave equation for the propagation of ¢ is 


n® " 5 


= 9° -+ ° 
A 2h? 


(4.4) 
C3 


Lamb’s (6) equation for the propagation of plane flexural waves in plates, 
obtained directly from the wave equation and the boundary conditions, is 


tanhrh =—s_ 4p*rs 


' Le2)2? 4.5 
tanh sh (p?-+s?)? (4.5) 


where r? = p?—n?/yc3, s* = p?—n?/c3, and r and s are either positive real 
or positive imaginary. It can be shown that (4.5) has an infinite number 
of solutions, each corresponding to a separate mode of vibration. The 
first mode, in which r and s are both real, reduces to Rayleigh waves for 
very short waves, and to the classical wave equation for very long waves. 
In the other modes s is imaginary, r taking real and imaginary values. 
Each mode above the first has a ‘cut-off’ frequency, given by hn > 4nc, 
in the second mode. For higher modes, the cut-off frequency depends 
on y; in the case y = 4, An > we, in the third mode, with the minimum 
values of hn considerably larger in the other high modes. The equation 
(4.5) is arrived at by considering only the components u and w of the 
displacement so that this is the equation with which (4.2) should be com- 
pared. The latter has two solutions for n (and therefore the phase velocity 
c = n/p) for each value of p. Corresponding values of p, c and the group 
velocity v have been computed for the first mode of (4.5), and for the 
mode of (4.2) in which c is smaller. The results are given in Tables | 
and 2, the ratio y being taken as 4 and 3 respectively, the corresponding 
values of Poisson’s ratio being 4 and }. The results show that provided 
ph < m, namely A > 2h, the plate thickness, the values of c and v com- 
puted from (4.2) and (4.5) differ only by one-tenth of one per cent at the 
most. An interesting result arising out of the computations is that the 
group velocity reaches a maximum of 1-01 lc, in the case y = 4 and 1-00lc, 
in the case y = 3, for wavelengths roughly twice the plate thickness, before 
falling to the velocity of Rayleigh waves for short wavelengths. 

Since in the first mode c < c,, and we are considering only ph < 7, 
we are dealing with frequencies such that nh < c,, so that only the first 
two modes of (4.5) are excited. In Table 3 there is a comparison between 
values of p, n, c and v computed from the second modes of (4.5) and (4.2), 
with y = 4. The agreement is fairly good, provided ph + 1. The cut-off 
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frequency for the second mode of (4.5) is given by hn > }ac,, whilst 
according to (4.2), 150/(y—1) 
—50y—10 ? 
From (4.6), An > 1-60c,, Fret en 8 the minimum value of hn by 
about 2 per cent for y = 4 and y = 3. 

If the expansion 


tanha = (1 ons 


(hn)? > Biy C3. (4.6) 





—_— — —_ 4.7 
3°15 315 (4.7) 
is used in (4.5) and only terms up to O(h*/A*) retained, an equation is 
obtained which is identical to (4.2), except for the last term which is 


- Be a 

807y*— 700 ry—140 4 
2100y(y—1)e$ 

The difference between this and the last term in (4.2) is 


h2n*/300(y—1)c4, 


a2 Qx* 172% 
oe 


which is small. 
The displacement ¢ is a rotation about an axis normal to the faces of 
the plate. Consequently (4.4) must be compared with the wave equation 
for the propagation of Love waves, which is (7) 
n2 (k+ y? 
= pty rae all 


Ce 


(4.8) 


for flexural vibrations, where k is a positive integer. If k > 1, hn > 3nc,/2, 
so that for the range of frequencies considered in this paper k = 0 is the 
only possibility. When k = 0, (4.4) and (4.8) are practically identical, a 
factor of 24 in the former being replaced by }7? in the latter. Thus the 
equation for ¢ in this theory yields a good approximation for the first 
mode of Love waves. 


5. Comparison with other approximate theories 

If all powers of (h/A) above the second are neglected then (3.16) reduces 
to 2 uh? 
via” Wy+piiy = 0. (5.1) 
This is the classical wave equation, and it can be seen from Tables | and 2 
that it gives good results only if the wavelength is longer than about 
twenty times the plate thickness. 

Mindlin’s (5) approximate equations are 


2 2 
[vs Fa aa” — Pop whey i 
3k*H 


h? 
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TABLE 1 
First mode 
ph c v Ca Va Cm em CMe Ome Ce Me 
o*100 0099 oO°197 0099 o'197 o°099 o°197 0°0g9 o°197 o°100 0°200 
0201 o19s5 0° 380 o'1gs 0° 380 orgs 0° 380 oO1gs 0 380 o°201 0402 
0°407 0° 366 0°665 0 366 0-665 0° 366 0664 0° 3604 06558 0°407 oS814 
0°620 0502 0°839 o*502 0°839 o°501 0°835 0°497 0824 
0°33g 0604 O934+, oho4 0°933 o*601 o927 o595 o’gIo 
o’g5i 0°644 o961 0644 0961 o-641 0°953 0633 0°934 
1*290 0°734 1002 0°734 1°OOI 0°729 o’ggI o-717 0-967 
1°748 0806 loll 0806 I‘O10 0°799 0°997 0784 "970 
1°g75 0°8350 1*009 0829 1-008 0-822 0994 0805 0966 
2°439 0863 1°000 o'862 1*000 0354 o"990 0835 0957 
2898 0884 0999 07883 ogg! 0°373 0976 0853 0945 
3357 0°395 o°gd2 0°595 0984 0°337 0*967 0°566 0940 
4501 o’gIg o°g65 o-918 0970 0906 0°956 0883 0°928 
g'O4I 0°932 0°936 0938 0953 0926 0°939 ogo! O'913 




















In this table, Poisson’s ratio is 4, c and v are the exact values of the phase and group 
velocities for wave number p (p 27/A, where A is the wavelength), and plate thickness 2h. 
€q and v, are the phase and group velocities obtained from the approximate theory of this 
paper, while M1 and M2 refer to results obtained by Mindlin’s theory, Mindlin’s constant k 
being taken as the velocity of Rayleigh waves, and 7/(12) respectively. c, and v, are 
velocities obtained by the classical theory. The velocity of rotational waves c, was taken 
as | for the purpose of all computations. 


TABLE 2 


First mode 








ph sd ° Ca Va “M1 M1 CMe °m2 Ce Ve 
0201 o1dss 0° 360 Or1d55 0° 360 o185 0° 360 o185 0359 OIG 0379 
408 "349 0°637 0° 349 0°637 0° 348 0°635 0° 345 0°632 0° 355 0°770 
0625 452 o813 0°482 o*813 0°480 0807 0°478 o802 
o*850 0°554 o’gi2 0°584 oori o-s81 0902 0578 0894 
1082 o°661 0963 0660 0962 0 655 0°950 o°651 o"940 
1°315 o°717 0°955 o*717 0987 o'710 0972 o"705 o’961 
"558 0’760 0995 o’759 0°997 o'751 0980 0745 0968 
1*799 o"792 1°OO1 0°792 1000 0782 0-982 775 0°969 
2041 o’817 0°999 o°816 0°999 0805 o*980 0795 0°966 
2-526 o551 o’9g1 oSs! ogg! 0838 o'972 0°830 0958 
3°01 0873 oo81 0872 8©6arg84 o859 8 ©=—_: 0" 0849 =—-0"'949 
3493 0887 o’g7!I 0°887 0976 0873 0956 0-863 0-942 
5°886 ogl4 0938 0919 0955 0°902 0°935 o890 839 0-092 
10°62 o’gig o"g20 0"932 0944 o’ol4 o-g25 o-go2 o’gl2 























The notation in Table 2 is the same as for Table 1, but Poisson’s ratio is }. 


In (5.2) and (5.3) w is independent of Z, u and v vary linearly with Z, 
H is a rotation similar to { and k is an undetermined constant, the value 
of which is fixed according to the needs of the problem. If we now put 
w = w'eXP2-") and [ = C’e"t-™), the wave equation which is obtained 
from (5.2) may be compared with (4.2) and (4.5) whilst (5.3) may be 
compared with (4.4) and (4.6). Mindlin used two different values of k, 
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one equal to the velocity of Rayleigh waves and for the other he put 
k® = w*/12. The first value of k was based on the fact that if this value 
of k is used the wave velocity of the first mode of (5.2) tends to the velocity 
of Rayleigh waves for high frequencies, as does the first mode of the exact 
theory. However, this leads to inaccuracies in the other modes, so that 








TABLE 3 
Second mode 
ph nh c v Ca: 1 Me Cy, Um Cys Unme 
0°244 1°636 669 O°512 6°83 0535 G@8gq 0546 671 O554 
0352 1700 «4°83 0689 4°93 9°722 4°98 0734 4°86 0747 
0439 1°766 4°02 0806 411 0848 415 0866 4°06 0876 
os16 183 = 3°55 BQO 3°63 0942 3°67 o9g62 3°59 o972 
0587 1897 3°23 O'955 3°31 «1016 3°35 =1°038 3°28 «17048 
0654 1°962 3°00 1°007 3°08 1°076 3:12 1100 3°06 «Il 
0779 27038 269 1°09! 2°77. «1°17 2°80 1°20 2°75 1°25 
1°O15 2°36 2°32 «1°16 2°41 «1°29 2°44 1°32 2°41 1°33 
1°359 277 2°04 =IQ 2°14 1°40 218 1°44 215 1°45 

















For this table Poisson’s ratio was taken as }, and n/2a is the frequency. The other 
symbols have the same meaning as in Tables | and 2. 


Mindlin used the value k? = 7?/12, which causes some inaccuracy in the 
first mode, but leads to the correct equation for the first mode of Love 
waves, and a correct cut-off frequency for the second mode of (5.2). 
Numerical results for Mindlin’s theory are compared with the other theories 
in Tables 1, 2, and 3, where ¢,,, and vy,, refer to Mindlin’s theory with k 
equal to the velocity of Rayleigh waves, whilst c,,. and v,. correspond to 
k®? = m*®/12. It can be seen from the tables that the theory developed in 
this paper yields slightly better results than Mindlin’s theory with k the 
velocity of Rayleigh waves. If k® = 7?/12, the phase velocity of the 
second mode is seen to be somewhat more accurate, at the cost of larger 
errors in the group velocity and also in the first mode. 

There is, however, little to choose between results obtained from the 
Mindlin equations and the equations of this paper. On the other hand 
Mindlin’s assumptions about the displacements are very wide of the mark, 
the discrepancies being hidden by a good choice of k. The fact that 
equations (3.16) and (3.18) compare so favourably with the exact equa- 
tions for low frequencies means that the choice of the shear stresses and 
transverse displacement made in this paper must be very close to reality, 
throwing more light on the behaviour of plates at low frequencies and 
also on the sort of accuracy to be expected from a Reissner type theory. 


6. Boundary conditions 
As with Reissner’s and Mindlin’s theories, three boundary conditions 
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are now required instead of only two demanded by the classical theory. 
It suffices to state the conditions necessary for various types of boundary; 
the proof is very similar to that of Mindlin or Reissner. 

It is assumed that the plate is bounded by a curve in the plane Z = 0, 
that the tangent and normal vectors at a point on the curve are s and n, 


and that n mykes an angle « with the axis of x. The conditions for a free 
edge are 


nz = ns = nn = 0. (6.1) 
If Q’ = ns—sn, 


I’ = ns+sn—2inn, 


pb’ = nz+i &, (6.2) 
then (3) O' = Q, 

IY == Te-™, 

os’ = pe-**, (6.3) 


From (6.1), (6.2), and (6.3) the conditions for a free edge are 

0+Te-*i« — 0, 
and re[pe-*™] = 0. (6.4) 
If the edge is horizontally clamped, 


w = D* = 0. (6.5) 
At a ‘simply supported’ edge we have 
w= u, = ne = 0. (6.6) 
With the use of (6.2) and (6.3), 
w = 204Te-%24[Pe2ia — re(D*e'*) = 0. (6.7) 


The initial conditions are that at t = t), w, D* and their time derivatives 
must be prescribed. 


7. Conclusion 


It has been shown that, provided hn < mc, the assumption that the 
shear stresses are parabolic and that the transverse displacement is of the 
form given in (3.2) leads to equations which give a good representation 
of the three possible modes of vibration. For higher frequencies there 
come into play other modes with nodal planes parallel to the faces of the 
plates, and no equation arrived at by the elementary plate theory can be 
expected to yield even a qualitative picture. 
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THE OSEEN FLOW PAST A CIRCULAR DISK 


By L. M. HOCKING (University College, London) 
[Received 20 November 1958] 


SUMMARY 


The streaming motion of a viscous liquid past a circular disk, fixed normal to the 
stream, is found on the basis of Oseen’s approximation by taking distributions of 
singular solutions of the equations of motion over the surface of the disk. The 
solution is found for low Reynolds numbers in the form of a series expansion in the 

teynolds number. The solution is also found for high Reynolds numbers, although 
the Oseen approximation is then invalid in the boundary layer and the wake, and 
the drag coefficient is found to be about twice the experimental value. 


1. Introduction 


THe Oseen method of solution for streaming motion of a viscous fluid 
past a fixed body is to suppose that the fluid momentum is convected 
with the free stream velocity instead of the local fluid velocity. This is 
an adequate approximation in the regions where the stream has been 
only slightly distorted by its passage past the body but it is seriously in 
error near the body and downstream of it, where the fluid may even be 
moving in the opposite direction to the stream. The solution of the Oseen 
equations can be found for all Reynolds numbers but only when the 
Reynolds number is small is the solution physically significant. It is of 
interest, however, to see what form the Oseen solution takes at large 
Reynolds number so the flow past a circular disk held normal to the 
stream is found here for both low and high Reynolds numbers. 

In a previous paper (1), a doublet technique for solving the Oseen 
equations for axisymmetric flow was developed, following the method 
used by Bairstow, Cave and Lang (2) for two-dimensional motion. This 
method consists essentially of finding a doublet solution of the equations 
of motion and regarding the complete flow past a body as that produced 
by a distribution of the doublets over the surface of the body, together 
with an irrotational flow. It is then possible to find an integral equation 
for the density of the distribution and hence to determine the flow pattern. 
In applying this method to the flow past a circular disk, difficulties arise 
from the coincidence of the two faces of the disk which suggest that 
doublets of higher order must be used. The way in which these doublets 
may arise from doublets of the first order can be seen by supposing the 
two faces of the disk to be separated by a small distance and the strength 
of the doublets at corresponding points on the two faces to be m-+-n and 

[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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m—n, say. Then as the two faces are brought together, a doublet of 
strength 2m is formed, unless n increases to infinity as the distance 
between the faces decreases to zero, in which case a doublet of higher 
order is formed in addition to the doublet of strength 2m. Distributions 
of the two kinds of doublets on the disk are used in what follows instead 
of the integral equation found in (1). 

Both cylindrical polar coordinates (x,r,«), and spherical polar coordi- 
nates (#, 6, «) are used in the following analysis of the flow past an 
axisymmetric body, the axis being given by r = 0 and by 6 = 0; « is the 
azimuthal angle. The Oseen equations for the flow of a stream of velocity 
U, density p, and kinematic viscosity v are satisfied if the axial and radial 
velocity components, u and v, and the pressure p, are given in terms of 
two functions ¢ and y by the relations 
mR. FS A 


“= Te abe ca wal 
eax 2k Ox 


xX 
l ex 
or | 2k er’ 


p = pl 


where 2kv = U; ¢ is a harmonic function and y satisfies the equation 
V2y = 2kéy/ex. Introducing a stream function defined by u = —éyb/r ér 
and v = éy/réx, the first-order doublet is defined by 

ys = (1+ cos 6@)[1—exp{k R(cos @— 1)}]/2k, (1) 
which consists of the stream function for a source and that corresponding 
to the simplest solution for x which is zerc at infinity, the magnitudes 
of the two parts being chosen so that there is no creation of fluid at the 
singularity. The second-order doublet can be found as the limit of two 
first-order doublets of opposite strengths as the distance between them 
tends to zero and their strength tends to infinity, as in the formation of an 
irrotational doublet from a source and a sink. The stream function for 
the second-order doublet is 


b = sin*[1—(1+kR)exp{kR(cos @—1)}]/ R. (2) 

If distributions of doublets of both kinds are used, their densities 
will not necessarily be uniquely determined for it was shown in (1) that 
the vorticity is completely determined if it is known over the surface, 
and various combinations of the two kind of doublets will give the same 
vorticity on the surface and hence everywhere. To the motion given by 
the doublet distributions. an irrotational motion may be added. But the 
doublets themselves contain irrotational parts, (i, containing the stream 
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function for a source and (2) that for an irrotational doublet. To make 
the doublet distributions unique the further condition is added that no 
irrotational motion is required in addition to that given by the doublets. 
The use of this condition can only be justified a posteriori but, using it, 
the flow can be found and as the boundary conditions determine the Oseen 
flow completely, this choice of distribution is justified. 

The disk is of radius a and lies in the plane 2 = 0, and the densities 
of the distributions at the point (0,6,«) on the disk are m’(b) and M’(b) 
for doublets of the first and second orders respectively. The stream func- 
tion for a ring of doublets can best be found from the total velocity 
componerts at any point produced by the doublets forming the ring. 
If this is done for the rings covering the whole surface of the disk, the 
complete stream function is found to be 

an a a 
i | er [1—exp{k(x—p,)}]m'(b)b dbda+ 


27 a 


+ [ | 2. . ©0801) F1 _ (14 kp,)exp(k(a—p,)}]M’(b)b dbda, 
a3 Pi 
where pj = x*+-r?+ 6?— 2rb cos a. 

The velocity components at the point (0,r,0) on the disk are 


27 a 
— a ke) m'(b)b dbda+- 
J¢ =p 
ules 27 a 


+{ j= L-+kp+kp*)exp(—ke) 915) dda, 
p? 





0 6 
27 a 
=f ’ (r—b cos x){1—(1+kp)exp(—kp)} m'(b)b dbda— 
kp? 


27 a 
[  (r—b cos a)(kp+-k*p*)exp(—kp) 


j M'(b)b dbdx, 
p 


~ « 
0 0 


where p* = r?+6?—2rbcosa. By changing the origin to the point where 
p vanishes and taking polar coordinates p and f, say, in the plane of the 
disk the integrals in (4) near p = 0 are proportional to ff (1/p)p dpdB and 
are finite. The second integral in (5) is proportional to {{ (cos 8/p)p dpdB, 
and this vanishes except when the point p = 0 is on the edge of the disk, 
as 8 can then only range between —}m and +42. To make v finite it is 
therefore necessary to add the condition M’(a) = 0. 
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The boundary conditions are u = —U, v= 0 on the disk. The 
expression (5) for v can be written 
P * —k 
al-l I 1~exp(~ke) 1(b)b dbda a+ | [=P P) yy’ (byb dba, 
, op 2kp 
and if the equations are made non-dimensional with a and U the units 
of length and velocity, the equations to determine the densities are 
- 3 K 
—i —4P) nib) dbdx+ 
2p 


/ @ 
0 0 





; ff 1—(1+ Kp+ K*p*)exp(— 


- Ke) wipys dbdx = —1, (6) 


~*~ « 
00 


j | Pr — =P) m\byb dbda— 


— Kp) M(b)b dbdx = constant, (7) 


(8) 
where m’ = Um/a, M’ = UMa, and ka = K, and 6 and p are now non- 
dimensional lengths. The Reynolds number Re is 2Ua/v = 4K. 

Goldstein (3) showed that the drag on a body could be expressed in 
terms of the flow across a large sphere produced by the function ¢ of the 
Oseen solution, i.e. the irrotational part. The contributions to ¢ made 
by the second order doublets produce no flow across a large sphere sur- 
rounding them so that the drag is 4mp(total strength of the first order 
doublets), as found in (1), » being the coefficient of viscosity. In terms of 
the density m(b) the drag is given by 


1 
drag = 8n°.Ua | m(b)b db. (9) 


2. The solution for low Reynolds numbers 


For small values of the Reynolds number, the equations (6) and (7) can 
be expanded in powers of K. If the densities m(b) and M(b) are expressed 


= m(b) = folb) + Kf,(b)+K*fq(b)+-... 
M(b) = go(b)+ Kg,(b)+K%9,(b)+..., 


the coefficients of successive powers of K in the expanded equations give 
a series of equations for the functions f, and g,. The first three equations 
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arising from (6) are 


27 1 

[ (0s dbdx = 2, (10) 
* Pp 

0 


0 
27 1 


27 
[ [a2 dbda— } [ so1o dbdx=9, (11) 
p J. 
0 


~ * 


0 0 0 


27 1 27 


on 1 
a [ | J:lb), dbda- [ Afi 0)—rfald) 20 dbda = 0 
* ¥ =p 4 bf <p 
0 0 


and the corresponding equations from (7) are 
27 1 


[ [ § f9(b)b dbda = 
00 
27 1 a 


[ { 27.0) af(0y}0 dbda— f [ 26 abn 
> > 4 * » Pp 
0 


0 0 0 
27 1 27 1 
rf ap ae ‘ 9 bd 9,(5) 
| A {p*fo(b) — 3pf(b) + 6f.(b) + 12g9(b)}b dbda— b dbdx 
J a ta «€*P (15) 
0 0 0 0 
where Cy, C,, and C, are constants. Equation (8) gives the conditions 
g,(1) = 0 (« = 0,1, 2....). (16) 
It is clear that the equations can be solved successively if equations of 
the type on 3 
TP ea 
fi ar dbdx = hi(r) 
J. P 
0 0 
can be solved. This equation is exactly that giving the potential h(r) on 
a circular plate on which the charge density is f(b). If A(r) is an even 
polynomial in r, it can be expanded in terms of the Legendre polynomials 
P,,,{,(1—r?)} and the value of f(b) when h(r) = P,,{.(1—r?)} is given by 
Bateman (4) as 


‘ > r > s __ pay) 
f(b) a “= o Wo ee Py iv b My (17) 


Thus (10) gives 


and this result substituted into (11) gives 


4 
b) = - aa 
A) 7/(1—6?) 


Equation (13) merely determines the constant C). 
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To ee the solution further, it is necessary to evaluate integrals of 


the type 1 f p"{(b)b dbdx (n > 1), for various known functions f(b), and 


66 


for this purpose the following result is needed. If 
27 1 


Fir) = | | Soyo dba, 


then {+ inl -f i( qf, S \y abda, (18) 
at rdr a * 


provided that in (18) the integrand is continuous except at a finite number 
of points and the integral is uniformly convergent. 

Under these conditions, the order of the integration and differentiation 
can be reversed so that 





a?F dF + ee ES d YY, 
dr? ' rdr- f i p dp _e 


0 0 





-f es ON adn feces fans £\ dbda, 


using (r—b cos a)® = p?—b?sin’a. The second integral is 


27 
* 


| eos (47) \0 doa = {ites “| bdo 
lr pdp' pdp}} Ji r pdelo 


and ss integrand is zero except possibly when 6 = r. If the integrand is 

finite there, the integrand has at most a removable discontinuity at this one 

point so that the integral exists and again is zero. This proves the result 

which will be used when f(p) = p”, and the above proof is valid for n > 1. 
In finding g,(b) from (14) it is necessary to evaluate 


27 1 
F(r) = | { pfo(b)b dbda, 
00 


where f,(b) has the value found above. Using the result just proved, 


id atop j jis fol), abda 


by (10), so that F(r) = 4774+ A+ Blogr. When r = 0, p = 6 and 


27 1 op? 
ron} | an 
00 
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so that A = 1, B= 0. The equation for g,(b) is then 
; { 1 } 
| | I) 5 dbda = —Jr?-+C", 
tv 

0 0 
where C; is a constant. Solving this equation in the way indicated above 
gives 


J 2— 3b? const. 
62,(1—62) * /(1—6?) 


Jo(5) 


and the unknown constant can be determined by the condition g,(1) = 0, 
giving g)(b) = ,(1—*)/27*. By using the methods described, the solution 
can be extended to any power of K. The drag coefficient can be found 
using (9) and, when the Reynolds number Re = 4K, is given by 

20-37 , 


“p Re {1-+0-1591 Re—0-0009 Re?—0-0003 Re*--...}. (19) 

The Stokes solution is given by the distribution m(b) = f,(b) of first- 
order doublets with no second-order doublets and gives the first term in 
(19). The next term was given by Oseen (5). If the values of Cp given 
by (19) are compared with the experimental values for falling disks 
obtained by Schmiedel (6) they are found to be in good agreement up to 
about Re = 5. 

The nature of the flow immediately downstream of the disk can be 
determined by finding the velocity along the axis. Including the terms 
of order K only, the velocity at a point distance x from the disk is found 
sty “u— u!?4 “(! +72)}(tan i: .... 4) UKs. 

a 7° 4/) z?+-1 a 

where u is now the actual fluid velocity and not that relative to the 
stream. When K 0, uw is positive for all x, but if K is not zero, the 
leading term in the expansion for small x is —}Kz, so that the velocity 
ciose to the disk is towards the disk. There is therefore always an eddy 
behind the disk, even for very small Reynolds numbers. The stagnation 
point through which the dividing streamline passes is at a distance equal 
to the radius of the disk downstream for Re about 2-5. It was shown in 
(1) that in the flow past a sphere to the same order of approximation 
there is no eddy behind the sphere; it is of course to be expected that an 
eddy will more readily form behind a disk than a sphere. 


3. The solution for high Reynolds numbers 


The Oseen approximation is valid, even at high Reynolds numbers, 
outside the boundary layer and the wake. Although quantitative results 
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derived from the solution of the Oseen equations at high Reynolds 
numbers are bound to be of little significance, it is of interest to find the 
nature of the solution to see if any light is thrown on the exact solution 
of the Navier-Stokes equation. This forms the subject of a paper by 
Stewartson (7) which includes the Oseen flow past a sphere at high 
Reynolds number. 

For large values of the Reynolds number, the exponential term in the 
stream function for the first-order doublet (1) is negligible except in a 
narrow paraboloidal wake behind the doublet. The flow produced is 
radially outwards except in this wake where there is a strong inward flow 
almost parallel to the axis so that there is no total flow of fluid from the 
doublet. The Oseen flow past a body at large Reynolds number will like- 
wise consist of an irrotational flow over all space except in the wake and 
in the boundary layer on the front surface of the body. 

In finding the flow past a circular disk at large Reynolds numbers, (6) 
and (7) must be solved for large K. It is supposed that K is sufficiently 
large for exp(— Ap) to be negligible when p > c, where c is small compared 
with 1. Since the exponential functions occurring in (6) and (7) are 
negligible except within a small circle surrounding the point p = 0, the 
integrals can be evaluated by changing to polar coordinates p and 8 with 
origin at a point distance r from the centre of the disk and by neglecting 
the variation of the density functions over the small circle. Thus 


27 1 


27 ¢ 
©xP(— KP) (by dda = m(r) [ [Po dpa 
P P 
0 0 


~ « 


0 0 
= 2nm(r){l—exp(—Kc)}/K = 2nm(r)/K 


since exp(— Kc) is negligible. The value of the integral is halved if r = 1, 
as the point can then only be surrounded by a semicircle, but this edge 
effect is ignored in what follows. Using (18) to transform the second 
integral in (6) and evaluating the integrals of the exponential functions, 
(6) and (7) become 

20 


1 
am(r), dj{d M(b) ps 
“rT +aA"s) {J ee 
0 


0 


27 1 
; | | melt dbda—2nK M(r) = constant, 
00 


where m/K has been neglected in comparison with m and similarly for M. 
Writing wm/(r)/K = f(z) and M(r) = F(x), where r? = 1—y?, the above 
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two equations and (8) become 


se ae d 1—p? F(a) ag — 20 
Swi waa rm in), J p 24 
| f(r) dS —4n*®F(u) = C, (21) 

p 
F(0) = 0, (22 


where y, is the value of » for any point of the disk, dS is an element of 
surface on the disk and p is the distance between two points given by 
and p, at an angular separation a. These equations show that f and F 
are of the same order in K. 

If we suppose that - 
fu) = ¥ Am Panl)/H, 


a form that is suggested by the solution for low Reynolds numbers, (17) 
shows that 


a dS = ¥ Ay, 7%, Pan(t), 


m=0 


ie oe 2m—1 
where x need aon 


ss 2.4.6....2m 
The value of F(y) is given by (21) as 


x 


«o 
F() > jA,, oe Py,,() i > }A,, es 
m=0 m=0 
the value of the constant being determined by (22). To evaluate the 
integral in (20) it is necessary to write F(x) in the form 


F(u) => B,, P; (He) be, 


so that equating the two forms of F(,), 


> By Poy (H) mit > 1A, a Got 1)F ‘2m +(e) + 2mP,,, _ (mu u) 


4m-+-1 


— > 44,08, Py). (23) 


m=0 


This equation has to be satisfied only over the half-range 0 < p < 1. 
The expansion of the odd function P,,, ,,(«) in terms of the even Legendre 
polynomials over this half-range is equivalent to finding the expansion 
of the even function /P,,,.,(u)|, |u| <1, in a series of Legendre poly- 
nomials, This can be done in the usual way using the orthogonality 
relations and, when substituted into (23), the coefficients of P,,,(4) on 
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both sides can be equated, giving} 


+ (4m-+1)a,, . 3 
— (—] tl a —j])" 
By ( ) 4 2,' ) a, A,X 


te Baden (2n+1)? = 
* |(2m—2n--1)(2m+2n+-2)(4n+1) 
(2n)? aS 1 
~ (2m—2n+1)(2m+2n)(4n+1)  (2m—1)(2m+2)f° 
The integral in (20) with the form of F(u) given above is a series in 
P,,, (uw) and 
a —p* “YP ee —_ 2m(2m+1){2M Pym (He) +-(2mM4+-1)Pom_1(4)}. 
wdp\ p dp} ™ (4m-+1)p* 


thus (20) becomes 








(24) 





2m(2m+1) 


y A,, »*P, + > wot B ——— 
6 mbt om f) m m 4m+ 1 


m 


{2m Pony iq(t.)+(2m+ VI) Pom—(H)} 


= §F(4)+8P(x), (25) 
after using the form of f(u) and multiplying by u*. The even function 
u?P,,,(u) must now be expressed in a series of odd Legendre polynomials 
over the half-range 0 < » <1. Forming the expansion as before and 
equating the coefficients of P,,,,,(u) in (25), another set of equations 
results, relating the coefficients A, and B,. If the values of B,, given by 
(24) are substituted in this set of equations, the equations to determine 
the A, are found to be 


m=0 


; 0-6 (m = 0) 
> Ant finn + (2M+ 12M) Gm,» +(2mM+2)(2M+-3) "Ym sa.n} = o (m = 1) 
0 (m> 1) 
where (26) 
finn = (—1)™*"(4m +-3)(2M+ Lot, op, X 
fy (2n+-1)? 
* | (2n—2m+1)(2n+ 2m +4)(4n+ 1)(4n+3) 
8n?24+ 4n—1 


~~ (2n—2m— 1)(2n+2m+2)(4n—1)ana3)* 


mts (2n)* ; 











(2n— 2m— 3)(2n+-2m)(4n— 1)(4n-+- 1) 
a (2n-+1)? r 
™\(2n—2m+1)(2n+ 2m-+2)(4n+ 1) 
(2n)? 
~ (2n—2m—1)(2n+ 2m)(4n+ 1) 
Ti 





Im,n rx i(— 1)™+" 1208, x 





1 
+ Genie FBP 
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The terms in the matrix of coefficients of A,, decrease quite rapidly away 
from the leading diagonal so that the solution can be found approximately 
by taking the first few equations. The drag coefficient in terms of the 


function f() is 1 


Cp = 8 | f(ulu du = 8Ag, 
0 

and if A, is found using the first 5 equations of the set (26), Cp = 3-6. 
The experimental value of the drag coefficient, for Reynolds numbers in 
the range where the fluid motion is not turbulent, is given by Schmiedel 
(6) as between 1-5 and 1-8 for Re about 100-200. It is to be expected 
that the Oseen solution will give too large a value for the drag coefficient, 
even when the drag is entirely pressure drag, as in the present example. 
For the velocity of convection in the Oseen approximation is the main 
stream velocity, which is everywhere greater than the true convection 
velocity. Other solutions of the Oseen equations at high Reynolds numbers 
have been given by Pierey and Winny (8) for the flat plate, using doublet 
solutions. This solution and that by Stewartson (7) for the sphere, both 
give drag coefficients which are considerably larger than the experimental 
ones. 

In describing the flow produced by the distributions of doublets on the 
disk at large Reynolds numbers, at points some distance from the disk, 
the doublets may be regarded as concentrated at the centre of the disk, 
so that the flow is produced by a first-order doublet of strength 2K A, 
and a second-order doublet of strength 27B,, the strengths being given 
by the total strengths of the surface distributions. The stream function 
for points not in the vicinity of the disk is 


y/Ua® = A,(1-+cos 6)| 1—exp{K R(cos 6—1)}]+ 
+ By sin?6{1—(1+ K R)exp{K R(cos @—1)}]/R, 
where A, = 0-449, By = 27 B, = 0-047. 


The flow outside the wake is given by the irrotational parts of this 
stream function together with the stream velocity. Inside the wake the 
exponential terms are not negligible and there is a vortex ring behind the 
disk. The velocity along the axis is U(1—2KA,/x+2K B,/zx?), so that, as 
K is large, the stagnation point is at a distance 2K A, along the axis, i.e. 
at a distance 0-22 Re times the radius of the disk downstream. As the 
Reynolds number increases, the stagnation point recedes to infinity and 
the wake narrows until it occupies only the cylindrical region behind the 
disk, of equal radius. 








THE OSEEN FLOW PAST A CIRCULAR DISK 475 


Acknowledgements 

My thanks are due to Professor T. V. Davies for his help and encourage- 
ment in the production of this paper. It is an abridged version of part of a 
Ph.D. thesis of the University of London. 


REFERENCES 
1. L. M. Hocxinc, Mathematika 5 (1958) 134. 
. L. Barrsrow, B. M. Cave, and E. D. Lane, Phil. Trans. Roy. Soc. A, 223 (1923) 
383. 
S. Gotpstern, Proc. Roy. Soc. A, 123 (1929) 216. 
H. BaTeMAN, Partial Differential Equations (Cambridge, 1944) 447. 
C. W. OsrEen, Hydrodynamik (Leipzig, 1927) 205. 
. J. ScHMIEDEL, Phys. Z. 29 (1928) 593. 
K. Stewartson, Phil. Mag. (8) 1 (1956) 345. 
N. A. V. Prercy and H. F. Wixny, Proc. Roy. Soc. A, 140 (1933) 543. 


nh 


DIM Ae 














NECESSARY CONDITIONS FOR 
OPTIMAL ROCKET TRAJECTORIES 


By D. F. LAWDEN (University of Canterbury, New Zealand) 
[Received 16 October 1958. Revise received 25 February 1959] 


SUMMARY 


The problem of transferring a rocket with minimum expenditure of propellant 
between two terminals in a given gravitational field and in the absence of aero- 
dynamic forces, is expressed as a problem of Mayer type in the calculus of varia- 
tions. The known solution to this problem is employed to yield a set of necessary 
conditions to be satisfied by the optimal trajectory. As an illustration of the general 
results, the particular case, when the field is uniform, is solved. 


1. Introduction 

Tue problem of determining rocket trajectories of minimum propellant 
expenditure in the absence of dissipative forces is of fundamental im- 
portance to the science of astronautics. In previous papers (1, 2), condi- 
tions necessarily satisfied by such a trajectory have been derived and 
then applied to the particular problem of deciding the optimal manceuvre 
by which a rocket vehicle might be transferred between circular orbits 
about two planets (3). Cicala (4) has remarked that all such questions 
relating to optimal rocket trajectories may be expressed mathematically 
as calculus of variations problems of the Mayer type. A very general 
problem of this type has been formulated and solved by Bliss (5) and his 
solution is accordingly available for the present purpose. Miele (6, 7) has 
made use of this technique to obtain equations from which trajectories 
of minimal propellant expenditure may be computed for rockets and 
rocket propelled aircraft moving in the vicinity of the Earth and subject 
to aerodynamic forces. In this paper, it will be shown that the method 
suggested by Cicala leads to the same criteria obtained previously for the 
case of a rocket moving in vacuo, provided the motor thrust is supposed 
unlimited in magnitude. However, more general criteria, which are 
applicable in the case when the motor thrust cannot exceed a certain 
magnitude, will also be found. 


2. Problems of Mayer type 
The general Mayer problem may be stated thus: 
2,(t) (¢ = 1, 2,...,”) are n unknown functions of an independent variable 
t which are to satisfy the m first-order differential equations 
$;(t,x;,,%;) = 90 (7 = 1,2,...,.m <n), (1) 
[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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over a range of values of ¢ including the interval t, <¢t <t,. The values 
taken by these functions at the end-points of this interval will be denoted 
by 


Xi(to) = Xo t(h) = ta, (2) 
and must satisfy the p equations 
Wyp(to, ty, Xiq%q) =O (& = 1,2,.... p < 2n+2). (3) 


If J is a given function of to, t,, 24, 2, i.e. 
J = J (ty, th, Tio» Xy); (4) 


it is required to choose the functions 2,;, subject to the constraints (1) and 
(3), and the end-points fy, t, so that J is minimized. 

To obtain a rigorous solution to this problem, it is necessary, of course, 
to state at the outset the analytical properties of the functions ¢,, etc., 
and to specify the set of functions x; which are admissible for the purpose 
of minimizing J. However, apart from the statement that the z, will always 
be supposed continuous together with their first derivatives, the %;, 
except at a finite number of points, these assumptions will not be made 
explicit, since the object of this paper is to apply known results rather 
than to describe the most general circumstances under which they are valid. 
The reader interested in this aspect should consult Bliss (5). 

Let A,, Ag,..., A,, denote Lagrange multipliers depending upon ¢ and to be 
determined subsequently. Let F be the function defined by the equation 


F = \y4;. (5) 


Here, as elsewhere in this paper, the usual summation convention is 
operative in respect of repeated lower-case literal subscripts. Let v,, v9,...,v 


Pp 
be constant multipliers, also to be determined later, and let H be defined 
by the ati 

by the equation Ht = J+od, (6) 


Then, the functions x; which minimize J necessarily satisfy the second 
order differential equations 


eF d (5 


dx, aaa, 


J=0 (= 1.2.40), (7) 


at all points where the #, are continuous. These are the Zuler equations. 
Also, at the end-points t = ty, t = t,, it is necessary that 


eH oH. 
at, t oz,” = 0, (8) 
CH oH 


a tig ta = 9% (9) 
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oH be (=) =f (10) 
0 


OXig \OX, 


m+ (EF — 0, (11) 
OX sy OX} 4 
the subscript ¢ ranging over the values 1, 2,..., n. 

Together with the m constraints (1), the n Euler equations determine 
the m+n functions z;, A; apart from 2(m+n) constants of integration. 
The (2n+2) equations (8)-(11), together with the p constraints (3) and 
the 2m equations obtained by setting ¢ = ty, t = t, in equations (1), 
determine the p constants v,, the end points fy, t,, and the 2(m-+-n) integra- 
tion constants. 

The derivatives z; may be discontinuous, but at each such discontinuity 
the Weierstrass-Erdmann corner conditions must be satisfied. These require 
that the following expressions shall be continuous at the discontinuity: 

oF oF 


F—z (12) 


= tan ? 
Ox; Ox; 


i ranging over the values 1, 2.,..., n. 
Finally, if the ¢; are not explicitly dependent upon ¢, the Euler equations 
possess a first integral, viz. 


» 4 
F—z,— 
OX; 


= constant. (13) 


3. Fundamental problem of astronautics 

This is the problem of calculating the trajectory joining two fixed 
points D and A in a given gravitational field, along which a rocket may 
be navigated with minimum consumption of propellant. The rocket 
velocity will be supposed prescribed at the point of departure D and at 
the point of arrival A, but the instants of departure and arrival will, in 
general, be open to choice. 

If M is the rocket mass, v is its velocity and c is its exhaust velocity, 
all at time t, the equation of motion will be taken in the form 


Mv = F+ Me, (14) 


F being the resultant external force acting on the vehicle, not including 
the motor thrust. Then, if Ox, x, 2x, is an inertial frame of reference con- 
sisting of three orthogonal axes and if f,(7,,2,,23,t) (i = 1,2,3) denote 
the three components along these axes of the gravitational attraction per 
unit mass at the point (z,,2,,2,) at time t, the equations of motion of a 
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rocket moving in the field can be obtained from equation (14) in the form 


$; 8 — 5p BS =0 


15 


$; = M+£, == 0 


where the v, are velocity components, the J; are the direction cosines of 
the thrust direction, f is the rate of propellant expenditure, andi = 1, 2, 3. 


The identity ¢, = 8@+8+R-1=0 (16) 


must also be added to the constraints (15). 

For any particular rocket, 8 will be positive and bounded and thus 
0 <8 < B. Following Miele (7), we shall make 8 a monotonic increasing 
function of a parameter a (of no physical significance), such that as « 
traverses the range (—00,0), B increases from 0 to B. The functional 
relationship will be chosen so that d8/da = 0 for sufficiently large positive 
or negative a. Thus, the vanishing of d8/da will imply that either 8 = 0 
or 8 = B, i.e. that the motors are either shut down or are operating at 
maximum thrust. During the manceuvre, « will be some function of t. 

The positions of D and A and the velocity of the rocket at these points 
being specified, the boundary conditions take the form 


p; = To—d;, = 0 
Pi.g = T—a,; = O 
Pig = Vig—D, = 0 
Pisg = Vy—A;, = 0 


(¢ = 1, 2,3) (17) 


where the d,; and a; are the coordinates of the points of departure and 
arrival respectively and D,, A; are the velocity components at these points. 
The instants of departure and arrival, viz. to, t; respectively, may or may 


not be specified. If they are, then additional boundary conditions exist, 
viz. v3 = tb —-T = 7 (18) 
ty = 4,—T, = 0 


All the constraints have now been stated. There are eleven functions 
of t subject to these constraints but otherwise arbitrary, viz. the z,, the 
v,, the l,, M, a. These will be chosen so as to minimize the characteristic 
velocity of the manceuvre, viz. 


V = clog(My/M,), (19) 
M,, M, being the rocket’s mass at D and A respectively. 
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The Lagrange function F takes the form 


F A(.— - Bl; fi +A; .9(%—0,) +A,(M+B)-+A 


The eleven Euler equations are accordingly 


1%, 


Clr; 


c 
X; A; l;, 
A, we A 


0 (- “A, 1 +s) -. 
."": on 
where i 1, 2, 3. 
From equations (21) and (22) it follows that the A; (i = 1, 2,3) satisfy 
the three equations ef, 
A, =A. (26) 
COX; 
These three quantities A; have the transformation properties of a vector 
relative to rectangular cartesian frames.t This vector is termed the 
prime r (2). 
From the equations (23) it follows that 
» 
A; 5 Agl; (i , 2, 3), 
i.e. that the primer is always parallel to the thrust direction (unless 8 = 0, 
when there is no thrust). 
From equations (25), it may be concluded that 
ae (28) 


MM 


In the second event, by consideration of equation (24), it follows that 


either dp/da=90 or A, 


7 oe . 
Ny = GAs 


or =, = — iz M. (29) 


7 
Integration now yields the result 
gia constant 
Sate 
+ The f; are known to be the components of a vector relative to such frames and it 
follows that @f;/0x;, are the components of a second rank tensor. 


(30) 


. 
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and hence, by the second alternative (28), 
A, 1; = constant. (31) 
But the A, are direction ratios for the thrust, so that 
I, = Aj/J(AZ+-A§+A}). (32) 
It now follows from equation (31) that 
4. (A?-+A3+A§) = constant, (33) 


i.e. that the primer is a vector of constant magnitude. 

By expanding the solution of equations (26) in a Taylor series, valid 
in a neighbourhood of a point ¢ = 7, it may be proved that these equations 
cannot possess an integral in the form of equation (33). The details will 
be found in (8). The first alternative (28) is accordingly inevitable, implying 
that 8 = 0 or B, ie. that the motors must either be shut down or must 
operate at maximum power. 

If 8 = 0, the rocket coasts under gravity alone, the components of the 
primer satisfy equations (26), A, = constant, Ag = 0, and the /; are in- 
determinate. 

If 8 = B, then as equation (27) is valid, the thrust and primer are 
aligned so that equation (32) is valid. Equation (27) may accordingly be 
written as 


B 
Ny = Sy VARY AE+ AS) (34) 


and then, from equation (24), 
. B 
— apy it A+B). (35) 


Since the optimal trajectory comprises arcs of two types, viz. (i) null- 
thrust arcs and (ii) maximum thrust ares, it is necessary to consider the 
Weierstrass-Erdmann corner conditions to be satisfied at the instants the 
motor is energized or closed down. It will be found that the conditions 
for the continuity of @F/é%, are that the A; (i = 1, 2,...,7) should be con- 
tinuous at such an instant. Since A, = —2,, ete. (equations (21)), this 
implies that the primer and its first derivative are continuous. There is 
no condition on Ay. The condition for the continuity of F—#,@F /éz, is 


that =A jp Blt fi) Anat MB 


should be continuous. Since the f; and v, are necessarily continuous and 
the A, have been proved continuous, this last condition only requires the 
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continuity of 


y PX ls —A, B. (36) 


Consider an instant at which the motors are brought into operation at 
maximum thrust. Immediately prior to this instant, 8 — 0. Hence, im- 
mediately afterwards, A 

u BX, 1,—A, B = 9, 
‘ 


A, = M 


A, 1; ve LA+A2). (37) 

However, if the motors are brought into operation at the instant of 
departure, corner conditions are not applicable and condition (37) need 
not be satisfied. 

Similarly, equation (37) must be true at any instant when the motors 
are closed down, provided this is not the instant of arrival. 

Over any portion of the optimal trajectory lying in a time-invariant 
gravitational field, the first integral (13) of the Euler equations is available. 
This takes the form 


a Pee Hi) Aji Vv; + A,B + A,(7+-24 l3—1) constant. (38) 


Making use of the constraints (15), (16) and of the Euler equations (22), 
this can be written 
: c 
Afi —Ayri al Vu Ail, —As} constant (39) 
Condition (37) being satisfied upon entry into or exit from a phase of 
maximum thrust, it follows that the expression 


e 


is continuous through the value zero at such instants; A,, A;, v;,, and f; are 
also continuous. The constant appearing in the right-hand member of 
equation (39) accordingly takes the same value over the whole of any 
portion of the optimal trajectory which lies in a time-invariant field. 
Consider now the constraints (17), (18) imposed upon the conditions 
at the end points. The function H is first formed thus: 
H = clog(M)y/M,)+v;,(xo—d,)+¥;..(%j,—@,;) +0%.6(Vig— D+ 


+44 .9(¥ —Aj)+¥43(to—TM) + a(t: —T)), (40) 


where v,, = 0 if the time of departure is not specified and v,, = 0 if the 
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time of arrival is not specified. Conditions (8)-(11) can now be written 


down, viz.: es 
Vig tM Lit Vise ot M, = 9, (41) 
0 


; : c : 
Ma Mista trie On — yp = 0, (42) 
MM, 


¥i—Ajss.0 ° Vise —Aio = (43) 
VisgtAjis, . Visgt Aja — (44) 
c c 


iil hie 45 
u, yt" — 


six conditions relating to the end values of the 1; being void. 

The conditions (41)-(44) serve only to determine the vy; and do not 
provide additional information concerning the optimal trajectory. If, 
however, the time of departure is to be chosen so that the propellant 
expenditure is minimized, then v,, = 0 and, by eliminating v,;, v;,,, between 
equations (41), (43), a non-trivial end condition upon the optimal trajec- 
tory is found. This is 


; a 
A, %; +Ay 3H -+ M = (46) 
and must be satisfied at the point of departure. Similarly, if the instant of 
arrival is not given, then condition (46) must be satisfied upon arrival. 
In view of equations (21), condition (46) can be written 


* ° c 
A; 0,—A; v; —_ we = (, (47) 


A, is constant along an arc of null-thrust and is determined by equation 
(35) over an are of maximum thrust. Equations (45) show that at the 
end-points A, must satisfy the condition 


= 7 (48) 


This completes the enumeration of the conditions necessarily satisfied 
by an optimal trajectory. The simpler form to which these conditions 
reduce in a case of prime importance will be investigated in the next 


section. 


4. The case of unbounded thrust 

In very many astronautical problems, although B is finite, the thrust 
duration necessary to effect the required velocity changes is small by 
comparison with the times spent coasting along the null-thrust arcs (a few 
minutes by comparison with hundreds of days). It will accordingly be 
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possible, with very great accuracy, to regard the thrust as impulsive and 
to neglect the motion during a phase of maximum thrust, i.e. the optimal 
trajectory will be taken to comprise null-thrust ares meeting at junction 
points at which the impulses are applied. 

On account of the relatively short duration of a thrust phase, the 
components of the primer and their first derivatives, as determined by 
equations (26), will remain very nearly constant over this period and we 
shall neglect any variation in these quantities over such a phase. But it 
has been proved that these quantities are all continuous upon entry into 
and exit from such a phase. It follows that they are continuous across 
a junction point from the end of one nuil-thrust are to the commencement 
of the next. This is the first set of simplified conditions. 

Suppose that a phase of maximum thrust commences at t = h and 
terminates at t=. Then, integrating equation (35) over the interval 
(h, k), it will be found that during this phase A, increases by 

, a 
| vel AB-LA2) dt 


h 


, 

, A,X, 

al® jos+re+anl—e ( es __ 
jar vases) —*° |) ep AapAg) 


h 


it, (49) 


where an integration by parts has been carried out to yield the final 
expression. But, since condition (37) is to be satisfied at each end of the 
interval (h,k) (provided neither is an end-point of the trajectory), it 
follows from equation (49) that 


k 
7 5b A; X; 
M V(X +3 +23) 


dt = 0. (50) 
h 
Neglecting the variation in the A;, A; during this phase, condition (50) 
becomes ke 
A, A; ’ dt 


- = 0. 51 
VORP RAR) | ne 
h 


The integral is clearly non-zero and thus it follows that 
A; X; = 0 
at a junction point. 

The above argument is invalidated if the point at which the impulse is 
applied is an end-point of the trajectory and condition (52) is not neces- 
sarily satisfied at such a point. 

If an impulsive thrust is applied at the point of departure, during this 
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short phase the variation of A, is governed by equation (35). Neglecting 
the variation in the components of the primer during this phase, equation 
(35) can be integrated over the duration of the thrust to yield 


MA, = VAL+AE+AD (G7) (53) 


But, at the commencement of the phase, A, is given by equation (48) and at 


its termination by equation (37). It follows that at the commencement 
of the phase 


VOR+AR+LA8) = 1. (54) 
Since the components of the primer are direction ratios for the thrust 
when the motors are operating, this implies that they are, in fact, direction 
cosines. 

If the motors are not operated at the point of departure, A, remains 
constant at the value given by equation (48) until the first junction point 
is reached, when condition (37) is applicable. This implies that equation 
(54) is valid at this junction and the primer components are direction 
cosines for the impulsive thrust there. At the termination of the first 
impulsive thrust, equation (54) is still valid since the primer components 
will not have altered appreciably since its inception. But condition (37) 
is to be satisfied when the motors are closed down. Hence A, = c/M at 
this instant and the above argument can be repeated to show that equa- 
tion (54) is true at the second junction point. Further repetition of this 
line of argument proves that the result (54) is true at all junction points 
and condition (48) is then automatically satisfied at the point of arrival. 

It follows from the first set of constraints (15) that 


; © PA.f f. 
A, 6; = pyri ahi = yet ade (55) 


since, on a null-thrust are 8 = 0 and at a junction point A; = /;. Condition 
(47) can therefore be written 

ASi—AyY%; = O (56) 
and is to be satisfied at the point of departure if ft, is not specified and at 
the point of arrival if t, is not given. 

Over null-thrust ares passing through a time invariant gravitational 
field, the first integral of the Euler equations (39) is available, and takes 
the form A, f;—A;, ¥; = constant, (57) 
the constant being the same on all ares. If the time of transit between 
D and A is not pre-determined, condition (56) shows that this constant 
is zero. 
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To summarize, if the thrust phases are effectively impulsive, the com- 
ponents of the primer must satisfy equations (26) and, at each junction 
point, (a) A,, 2 
(b) A, A; 0, 


(c) A, must be direction cosines for the thrust. 


must be continuous, \ 


(58) 


These conditions imply that, at a junction point, the primer and its 
derivatives are continuous and that these vectors are orthogonal. 

If an impulsive thrust is applied at an end-point, condition (4) is not 
effective at that point. 

In addition, if the times of departure and/or arrival are not pre-deter- 
mined, condition (56) must be satisfied at the appropriate end-point. 

If the gravitational field is time invariant, equation (57) is valid over 
all null-thrust ares. Condition (6) then follows from the circumstance 
that the right-hand member of this equation assumes the same constant 
value over all such ares. For, A,, A;, f; being continuous at a junction 
point, it follows that, if »>, vj are the velocity components immediately 
prior and immediately subsequent to the impulsive thrust applied at such 
a point, 


But wi —o; oc l Asi 


Hence, condition (5) follows. 


5. Optimal manceuvres in a uniform field 

As a simple example of the application of conditions (58), consider the 
problem of transferring a rocket between two points D and A in a uniform 
gravitational field, the terminal velocities and the time of transit being 
specified. 

Since the field is uniform, the right-hand members of equations (26) 
vanish and the primer components are all linear functions of the time. 
Condition 58 (a) requires that the numerical coefficients occurring in these 
functions should not alter throughout the motion. Thus, if the primer 
is taken to be the position vector of a point P relative to an origin O, the 
locus of P will be a single straight line L. At a junction point, which is 
not the point of departure or the point of arrival, conditions 58 (6), (c) re- 
quire that OP should be of unit length and perpendicular to L. It follows 
that there can be no more than one such junction point. But, if such a 
junction point exists, the points of departure and arrival cannot be junc- 
tion points for, if they were, by condition 58(c) the corresponding posi- 
tions of P would have also to be at unit distance from O, which is im- 
possible. However, in general, it is easy to show that the manceuvre 
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cannot be accomplished by the application of a single impulsive thrust. 
It follows that the impulsive thrusts must be applied at D and A. 

The manceuvre is now uniquely determined by the boundary condi- 
tions. Only one parabolic trajectory exists along which the rocket can 
coast from )) to A in the given time of transit. Hence the impuise at D 
must convert the rocket velocity into that required for entry into this 
trajectory and the impulse at A must convert the velocity of arrival at 
A into that required at this point. 

Conditions 58 (a) and (c) have to be satisfied at D and A. The impulses 
at D and A being known, condition 58 (c) determines the primer vector 
at these points and hence the positions of the point P. The straight line 
joining these positions is then the line L. 

The conclusions of this section can also be obtained by employing 
elementary techniques as explained in (9). 
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SUMMARY 
For cases when the functions a(t) and y(t) have rational Laplace transforms, 
x x 
simple methods are presented for the evaluation of the integrals [ xy dt and f x* dt 
0 0 


in terms of the transform coefficients. The results are extended to the moment 
x 


integrals { try dt and f t'x? dt, where r iar 
0 0 


1. Introduction 


IN numerous design problems a fundamental part is played by integrals 
of the forms 


x x 
- 


l | ay dt, J ( x? dt. 

0 0 

For example, if in a control system of any type a(t) represents the instan- 
taneous error, i.e. the difference between the actual and the desired output, 
the integral J gives a ‘quadratic measure of error’ which is frequently used 
as a design criterion, and has to be minimized for optimum design. Again, 
if a(t) represents the current flowing in any (simple or complex) element 
of an electric network in response to a stimulus applied at time ¢ = 0, and 
y(t) represents the potential difference across the element, the integral / 
gives the total power dissipated in the element. Similar instances occur 
in other contexts. 

Since for linear systems design is normally carried out in the domain of 
the complex frequency variable p, the functions usually met in the first 
place are not 2(t) and y(t), but their Laplace transforms #(p) and 9(p), 
and in the majority of cases these are rational functions of p. The direct 
evaluation of the integrals may then be a very laborious process, involving 
(a) finding the zeros of the polynomial denominators of # and #, (b) obtain- 
ing x and y, the inverses of ¢ and #, by partial-fraction expansion or some 
equivalent process, (c) finding the product 2? or xy, and (d) integrating 
each of the possibly numerous terms in this product. 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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It is well known that some of this work may be avoided by using 
Parseval’s formula to replace the integral by one taken along the imaginary 
p-axis, the integrand being a product of transforms. The new integral 
can be evaluated by the method of residues, after the denominator zeros 
have been found, as in (a) above; but it was shown by Phillips (1) that 
the Parseval integral for J can be expressed in terms of the coefficients 
of the denominator, thus avoiding the often heavy labour of finding the 
zeros. Phillips did not give a closed expression for the value of the general 
integral, but only a rather complicated process for obtaining it. Later, 
Oldenbourg and Sartorius (2) gave a determinant formula for the integral, 
in a somewhat clumsy form and without any indication of the method 
of derivation. Recently, Babister (3) obtained similar formulae for J 
taking the differential equation of the system as starting-point. 

In the present paper a new and simple method of evaluation of both 
forms of integral is presented.t They are shown to be expressible as one 
of the unknowns in a set of simultaneous equations whose matrix involves 
the coefficients in # and 7. Thus the integrals can either be expressed 
formally in terms of determinants, or evaluated numerically by any of 
the usual methods for the solution of simultaneous equations. The labour 
of solution can be reduced by arranging the work so as to exploit the 
special form of the matrix, and this is discussed here in connexion with 
the Gauss—Doolittle method. A special method for J based on Routh’s 
algorithm is also described. A novel method for J and J which involves 
the use of continued fractions and requires less computation than other 
methods will be presented in Part II of this paper. 

Since for some control-system applications the ‘time-weighted’ measures 
of error p! 


A= [ ttdt, J,=[ @xrde 
0 0 
are advocated (4, 5, 6) in preference to the simple quadratic measure, 


both these integrals and the more general ones 


1 =j tay dt, I= | toy 


are also treated in the present paper, by simple extensions of the methods 
for J and J, and a further extension to deal with the higher moment 
integrals e 7 
I= | txydt, J,= { trx? dt 
* . . 0 a 
is indicated. 

+ The method used by Phillips for J does not seem capable of being generalized to deal 
with J. As far as the writer is aware, J has not previously been tackled at all. 


5092 48 Kk 
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2. Evaluation of J 
Let 2(t) and y(t) be two functions possessing Laplace transforms: 


E(p) ( ze dt, ifrep >a, 
0 


i(p) = | ye-'dt, ifrep > B, 
0 
where « and £ have their (algebraically) lowest possible values (i.e. they 
are the abscissae of convergencet of # and #7). Thent #(p) is regular in 
the half-plane rep > «, and g in the half-plane rep > 8. Suppose now 


that - LB - 0. (1) 


Then the singularities of #(p), which lie in rep < a, are all te the left of 
those of 9(—p), which lie in rep > —8 (and likewise the singularities of 
j(p) are all to the left of those of #(—p)). Now let y be any number 
between a and —f: a < y < —f, and let C be the line rep = y. Then 
C is to the right of the singularities of #(p), and provided #(p) tends to 0 
uniformly on the infinite semicircle to the right of C, and is absolutely 
integrable on C, C may be taken as a Bromwich contour in the inversion 


formula for x, viz. 





l a7 
x(t) ' E(pemdp (t > 0). (2) 
2m 
a 
At the same time, the transform 
( ye dt = 7(—p) 
0 
exists for p on C’, where re(—p) = —y > BB. Thus§ 
. -_ 
I y dt E(p)e” dp 
q 27 . 
0 Cc 
ee > 
¥(p) dp ye” dt 
2nt . . 
Cc 0 
oe Fe, ‘ 
=~ | =(p)i(—p) dp, (3) 
2m | 
5 


+ G. Doetsch, Theorie und Anwendung der Laplace-Transformation (Dover Publications, 
1943), p. 17. 
t Ibid., p. 43. 
§ The interchange of integration order can be justified under the assumptions made, 
x 
provided that f y\ev' dt converges, i.e. —y is greater than the abscissa of absolute con- 
0 
vergence of 9(p). 
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which is the form of Parseval’s formula needed here. Obviously x and y 
may be interchanged on the right-hand side. When Z and @ are rational 
functions, vanishing at infinity, with the poles of Z(p) all to the left of 
those of 7(—p), the formula is certainly valid. Since, then, 


F(p)j(—p) = O(p-*), 


the contour C in (2) may be completed by an infinite semicircle on either 
side. Taking it, say, to the left, and writing 


=p) = K(p)if(p),  H—p) = k(p)/9(p), (4) 
we obtain, by Cauchy’s residue theorem, 
I= > (Nif'9) pn, (5) 


where .V denotes hk, and p,; are the zeros of f(p), which at first we assume 
to be all simple. None are zeros of g(p), since by hypothesis all zeros of f 
(i.e. poles of Z(p)) are to the left of the zeros of g (poles of #(—p)). Thus 
f and g have no common factor. 

In order to evaluate the sum in (5), we consider the identity 


N uv 
 — +o, 6 
fg f g (*) 


Since g and & have no common factor, and the degree of NV, or hk, is less 
than that of fy, uf is the (unique) sum of those partial-fraction terms of 
N fg which belong to f; and similarly for v/g. Thus 


u SAL», 
f p 3 P—Pi 


It follows that u/f is, in general, of degree —1, and that J is the coefficient 
of | pin uf, if this is expanded in descending powers of p. Thus if 
f= fothp+--+hnP™ (Sm # %), 
J = Jor MU P+--+GnP™ Gn F%), 
hk = N= Not+Np+...+N,p” (N, # 0),F 
U = Ug t+ Uy, p+...+U,,-) p™-*, 


V = Ug $O Pt. + Uy P™, 


then I= u,_; /f,.- (7) 
Now (6) is equivalent to the polynomial equation 
fo+gu=N. (8) 


By equating terms in (8) we obtain a set of simultaneous equations for the 


+ In the present instance the degree y is at most m+n — 2, but all subsequent work applies 
equally well if N has degree m+n—1. 
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coefficients in v and u, which may be written 


Rw = N, 


where w is the column 


N is the (m+n)-rowed column 
(No. Mi, 
and R is the (m+-n)-rowed square matrix 


fo =O 


Sin 
In 
a . 2 “ed 
Thus formally we have = A’/f,, A (10) 
where A = det R and A’ is the determinant + daniel from A on replacing 
its last column by N. It will be seen that A is simply the resultant of f 








and g, and does not vanish since these polynomials have no common 
factor, by hypothesis. 

Finally, it is clear by continuity considerations that (7) and (10) will 
remain valid even if the zeros of f(p) are not all simple. 

It may be noted that the whole procedure above is essentially an 
algebraic process for evaluating the sum (5), and holds for arbitrary 
numerator NV, and arbitrary denominator factors f and g having no common 
factor, provided the degree of N is less than that of fg. In particular, 
the formula (10) will correctly represent the sum (5) even if the zeros of 
f are not all to the left of those of g. It will not, however, represent ‘the 
integral J, which in fact in such a case would not converge. For example, 
suppose x cos wt, y eft. 

Then i(p) = p\(p?+e), rep>O («a= 0), 
i(p) = 1/(p—B), rep > B. 


If 8 > 0, the integral [ xy dt does not converge. 
0 
If 8 < 0, ie. (1) holds, the integral converges to —8/(w?+§*). Also the 
Parseval integral (3) an 
—pdp 
2+ «y?)/ 2)(p+B) 


and the sum (5) becomes 


> ea Bo. So 
2p(p+8)  (B*+w?) 


p=+tiw 
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which of course holds for any value of 8, positive or negative. Again, we 
have ' 8 


io Pl 
A=|0 1 B|=o*+*, 
1 
0 | 
A’ =| 1 —1|=-—8, 
y ace 
and since f,, = 1, (10) gives the correct value, —f/(w?+-8*), of the sum 
(5) for all values of B. 
As has already been remarked, the roles of x and y may always be inter- 
changed. If this is done in the example just treated, we obtain 


a —pdp —p —B 


I — = = = > 
2m J (p—BYp*+w*)  (p*t+w*)pag (B+!) 
and with f = p—8, g = p*+w?, N = —p, we have 


|—Bp . 0 
A’=| 1 —B -1|=-8, 
paymnye § 0| 





giving J as before. 

In numerical examples, the equations (9) may be solved by any of the 
standard methods. Certain simplifications occur on account of the special 
nature of the matrix R. We shall illustrate the procedure by finding /, 
given that pP+p—3p+2 _ Ap) 
pap 8p pel Jip) 
EE os ak 8 oe ES | 
' p’+3p*+ 15p3+-36p*+40p+27 = g(—p) 


t(p) = 


The denominators of #(p) and 9(p) are both Hurwitz polynomials, so that 
the condition (1) is satisfied. Evaluating N(p) = h(p)k(p), we find that 
N = (—2, —3, 14, —11, 5, —8, 5, 4, 0). 
In the table which follows we use the Gauss—Doolittle method for solving 
(9). (The working in the Choleski method is very similar.) 
The last row yields the value u, = — 21/88, whence by (7), J = —21/88. 
In the general case the first n rows of multipliers consist of —f,/fo,..., 
J nfo arranged in a diagonal pattern, and the first n entries in the modi- 
fied (lower) uy column are repeated in the columns %,..., %,,_, in a similar 
diagonal pattern, which is, however, cut off horizontally after the (n+ 1)th 
row. The modified columns headed v,..., v,,_, contain merely the elements 
f, in the main diagonal. It is clear that in practice the columns %p,..., v,,_, 
may be omitted altogether. 
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My Ug Us 





Multipliers 





—-3 —-2 -! 
I -3 -2 —- 
—3 3 

=z 


I 
-2 
3 
I 

















3. Evaluation of J 

This is merely the special case of the previous integral with y = 2, but 
is best treated independently. 

First, the condition (1) now becomes « < 0, so that the poles of #(p) 
must lie in the (open) left half-plane. Thus, if #(p) = h(p)/f(p), f(p) must 
be a Hurwitz polynomial. 

Equation (3) becomes 


. 


l ao 
J=— | ¥(p)z(—p) dp, (11) 
2ri 


C 
where C is the path rep = y, a < y < —a (and may, for example, be 
taken as the imaginary axis), and in place of (5) we obtain 


J = ¥ (h(pyh(—P)/f'(P)f(—PN pve 


Clearly, if u/f is the sum of partial-fraction terms in h(p)h(—p)/f(p)f(—p) 
belonging to f(p), then the sum corresponding to f(—p) is u(—p)/f(—p), 
since the total sum is even in p. Thus in place of (8) we obtain 
f(p)u(—p)+f(—p)u(p) = h(p)h(—p) = 22, say, (12) 
and again J is the coefficient of 1/p in u(p)/f(p). 
If now we write 

f(p) = fot fipt--+SmP™ (Sm # %), 

u(p) = Upt+u, p+...+u,,-)p"", 

Lp) = Lyt+ Ly p?+...+ Lem. p?"—*, 
J is again equal to u,,_,/f,,- 


If the Gauss—Doolittle method is to be used for firding 1,,_,, we must, 
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in order to benefit fully from the special nature of the equations involved, 
arrange them thus: 





—_ Lem-s 


Once again, in practice, only the second group of columns on the left need 
be retained. The first [4(m+1)] rows of multipliers consist of —f, /fy, 

arranged in a diagonal pattern, and the first [4(m—1)] entries 
(excluding the initial zero entry) in the modified u, column are repeated 
in a similar diagonal pattern, which is however cut off horizontally after 
the next row. 

If m is even, the last column retained is headed w,,_,, and so yields J 
directly. If however m is odd, the last column retained is headed w,,_», 
and in this case w,,_,, and so J, are then found by using the last equation 
of the set, viz. i a 2 as (14) 

To illustrate, we calculate J when 
___ pP®—p>+3p*—2p?—Sp*+p+3  __ A(p) 
p+ p+ 6p + 5p*+ 9p + 6p?+3p+1 — f(p) 

We find L = }(9—31p?+ 47p*— 26p*— 5p + 5p!°+ p). 





«(p) 


uy Us Us 





—3 
—6 


nae —s/2 


— 5/2 
Multipliers 1/2 





—s 
—5 = 
—6 —§ 


s 256 
_ | 





6806 
ssaet 
454 
257 
82 














We thus obtain wu, = 257/2 and, by (14), u, = 129, whence J = 129. 
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Reverting to the equations (13), we obtain, after rearranging columns, 
the determinantal solution 


J = (—1)"*'!D/f,, D (15) 
fo 
fe fh fo 0 

D = hh fh h fo 


fom-2 . : . Sart 
in which f, = 0 (r > m), and D is obtained from D by replacing the last 
column by the column (Lo, Lg,..., Le,-2). D is of course the principal 
Hurwitz determinant for f(p), as well as the resultant of f(p) and f(—>p). 
It is certainly non-zero, since its vanishing would imply the existence of 
zeros of f(p) on the imaginary axis or in the right half-plane. 

We now present a special method for J, based on (15), which is of some 
interest, being an adaptation of Routh’s algorithm for determining the 
number of zeros of a given polynomial in the left and right half-planes. 
As is well known, this algorithm corresponds to a triangularization of the 
Hurwitz determinant D by a series of elementary operations. A modifica- 
tion of the algorithm leads to a similar triangularization of the determinant 
D in (15), from which J can be written down at once. 

To show how this may be achieved, we write 

f= (Foot Fon P+ Fos P*+ --)+ (Fir p+ Fis P+ --.). 
Then D == det M,, where M, is the (mx m) matrix 
Foo Foe 
F, —_— 
Foo Foz Fos 
Fis 


where 





L F. 


qm-1 
Here, F., = 0 (¢ > m), and q = 0 (m odd) or 1 (m even). 





Suppose now, starting with M,, we derive a sequence of (mm) ma- 
trices M,, Msg,..., M,,_,, by the following successive sets of modifications: 
for M,, subtract (Fpo/F,,) x rows 2, 4, 6,... of M, from rows 3, 5, 7,...; 


for M,, subtract (F,,/ Fy.) x rows 3, 5,... of M, from rows 4, 6 


for M,,_,, subtract (F),3 »—-3/ Fm—2m-—2) x row m—1 of M,,_, from row m. 


Here, Fy», Fy3,... are the leading elements of the newly-produced rows 
of M,, Mg,... respectively. 
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It is easy to see that, on account of the special form of M,, elements 
below the main diagonal progressively disappear in the sequence M,, and 
that the final matrix M,,_, is triangular, and of the form 
[Foo Fon Fos 
Fy As 

Fry 
Fy 
O 





F-1m-11 
Since the matrix determinants are invariant under the operations (16), 
we have 





D = F,, F,,, ... F-s.0-1° (18) 
Now M,,,_, is identical with the matrix formed by Routh’s algorithm for f, 
and may be obtained directly as follows. Let the row-vectors of M,,_, be 
denoted by ry, Fj,.--, Ty,—-;, and let S denote the operation of shifting the 
elements of a row one place to the right. Then starting with r, and r,, 
obtained from /, we derive ry, r,... in succession by the recurrence formulae 


Via — Sr;_,—B;¥; (a — 1, 2,...,m—2), (19) 


where 8; = F,_,,-,/F,;, F,; being the (i+1)th element of r; (i.e. the 
element on the main diagonal of M,,,_,). 


We now consider D. When transposed, its matrix M, ditfers from M, 
only in having the vector 


Vo => (Lo, Dog,.05 Lm-2) 


as its last row. It follows that if the operations (16) are applied to M,, 
except as regards its last row, the resulting matrix will be triangular except 
for its last row. If however simultaneously with the sets of operations (16) 
we operate on the last row in such a way as to make successive elements 
vanish, we shall obtain a triangular matrix. It is easy to see that this is 
accomplished by replacing y, successively by y,, Yo,..-, ¥,-1. Where 
Yiur = Ye —aF, (¢ = O,1,...,m—2). (20) 
Here, a; = H,/F,,, H, being the (+ 1)th, and in general the first non-zero, 
element of y;. The determinant is again invariant under these operations, 
so that D = Fey Fry... F-am-2 He-t- (21) 
We then have, by (15), (18), and (21), 
J = (—1)" on -a/fa- (22) 
For practical purposes the two sequences (19), (20) may be combined, 
a convenient ordering of row-vectors being 


Yo: Tor Yas Tis Yas---> Ym—a» Tm 
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where yp, fy, and r, are written down initially. In this scheme, each y-row 
is formed from the two immediately preceding rows of the scheme, each 
r-row from the two previous r-rows. 

We illustrate the method with the previous example. 

Yo 

To 
vi 
ry 
Ya 
tT, 
Ys 
Ts 
Ys 
T% 
Ys 
rs I 
Ys 43/14 
Tr, 1/42 
The multipliers «; and 8; have been written in line with the rows r; which 
they multiply. Using (22), we find J = 129, as before. 

It will be observed that a number of the entries in this scheme, at the 
ends of the rows y;, are repeated one or more times and, in contrast with 
the situation in the Gauss-Doolittle method, such repeated items need 
not be re-entered. Taking this into account, we find that in the general 
case the algorithm involves about 4m?+ 4m entries, 4m?+-4m multiplica- 
tions, 2m divisions. The corresponding numbers for the Gauss—Doolittle 
method are about {m?+-m, jm3+-}m?, bm?— 4m. The algorithm is obviously 
superior for, say, m > 5. A glance at the explicit formulae given by 
Phillips (loc. cit. (1), p. 370) for the case m = 7 will immediately demon- 
strate the comparative simplicity of either of these methods. 


4. Evaluation of J, and J, 


For some purposes, as was mentioned in the introduction, the time- 
foo] 


weighted measure of error, J, = | tx? dt, is preferred to the unweighted 
0 


measure J. The adaptation of the foregoing methods to the evaluation of 
this integral and of the more general integral J, will now be discussed. 
As is well known, the formula (3) may be generalized into 


© 
fa 


e-@ xy dt = = [ E(p)i(q—p) dp, (23) 


0 C 


where C is the line rep = y, anda < y < req—f. Assuming, as in section 
2,a+B8 < 0, let y be chosen in the range « < y < —f. Then q may be 
situated anywhere in the half-plane reg > y+, which includes q = 0, 
since y+B < 0 
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The left-hand side of (23) is the Laplace transform of xy with respect 
to q. Since it exists and is regular in a region including the origin in the 
q-plane, we may expand both sides of (23) in ascending powers of gq and 
equate coefficients. This gives 


a . —l “er 
I, = coefficient of g in = | Hova—m dp = —. ani | HP) p) dp, 
(24) 
r, using (4) again, = a = (i, ve) dp, (25) 


each term of which may be evaluated by the method of section 2. How- 
ever, the second term here may be of rather high order, and an alternative 
method avoiding the increase of order, though of course involving some 
increase of computational work, will be described. 
By (24), we may write, retaining only terms of degrees 0 and 1 in q, 
ee St | Fama qe) 4 
fig—99') 


h(k+-qk’ ) 
fig+q’) 


Now suppose that polynomials u, io, v, v» can be found such that, as 
far as terms of order q, 

f(v+qv?)+(g+¢g')(utqu®) = h(k+qk’). (27) 
It will then follow as in (7) that 


I, = «_, If... (28) 

To determine w‘!)_,, we note that (27) is equivalent to the two equations 

fo+gu = hk, = N, (29) 

fv™+gu® = hk'—g'u, = N®, say. (30) 

Here, (29) is identical with (8), while (30) is of the same form, so that J, 

can be determined from (30) after finding « from (29) and inserting it in 

the right-hand side of (30).t 

If the Gauss—Doolittle method is used, the additional work in finding 

I, beyond that involved in finding J is not great. For example, to find J, 
for the case treated in section 2, we have, by back-substitution, 

us = —21/88, wu, = 1/88, u, = —21/44, 


+ Equations (28)-(30) can in fact be obtained directly from (25). For if u and v satisfy 
(29), and N“) is defined as in (30), the integrand in (25) is N“/fg +vg’/g*. The integral of the 
second term vanishes since its poles are to the right of C, and the result (28) follows at 
once by (7). However, the method in the text has the advantage over this method that it 
can be readily extended so as to deal with /,, J,,... a8 will be seen in section 5. 


= coeff. q in | Fass 


(26) 
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so that 

N® — (—928, 888, 1294, —394, — 1421, 505, 1665, — 105, 0)/88 
Thus, using the same multipliers as before, the modified last column 
becomes 


(—928, 1816, 2262, —6248, — 4663, 17572, 226933, 


1, — 194635, —747})/88, 


3736/5. 88 
- 3736/5 . 88 ~ 467/968. 
88/5 


We turn now to the integral J,. We have a < 0, and 


giving J, = uf) (by (7)) = 


J, = coeff. of q in — | E(p)%(q—p) dp, (3i) 


- a 


where C is the line rep = y, a < y< x, and reg > y+a. This may 
of course be evaluated in the same way as /,. If we use the notation 


f*(p) orf*=f(—p), f'*=f'(—p) = —(f*y, ete. 
(29) and (30) become fu+ftu = hh*, (32) 
fvP+f*tu® = hh*’—f*'u. (33) 


In (32), the right-hand side is even, and v must be u*. Thus wu can be 
found by the Gauss—Doolittle method as in section 3, using back-substitu- 
tion. (Note that the modified columns omitted, under up, wg,..., contain 
only the elements f, in the main diagonal.) However, on turning to (33) 
we find that the right-hand side is no longer even, so that the same simple 
method no longer suffices. This situation, which arises because the inte- 
grand in (31) is not even in p, can be remedied by the simple device of 
changing the variable of integration. If we replace p by p+ 4q, the integral 
becomes 


| &(4q+p)z(4q—p) dp, 


Cy 


where C; is the line rep = y—re }q, and the integrand is now even in p. 
Thus we have 
. a 
J, = (—4) coeff. of g in Mp gn"(P—) q) dp. (34) 
2mt J flp+asf*(p—9) 
N ow if q lies in the region y+a < reg < 2(y—a) (which includes g = 0), 
C, can serve as the path C in (11). It follows, as for J,, that if, up to terms 
of order gq, 


(f+af’)(u*+-qu*)+(f*+qf'*)(u+qu®) = (h+qh’)(h*+-qh’*), 
fu*+f*u = hh* = 2L, (35) 
fu®* +f*u) = hh’'*+h’h*—f'u*—f'*u, = 21, (36) 
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where L and L® are both even, then 
Jy = — Un 1/2fmi (37) 
u‘)_, is found from (36) after uw has been found from (35) and inserted into 
the right-hand side of (36), which, we note, gives 
L® = even part of (h’h*—f'u*). (38) 


As a simple example, we shall evaluate ( tx? dt, given that 
0 


E(p) = (12p3+- 16p?+ 10p+3)/(p+1)(2p+ 1)(2p?+ p+1). 


Here, 
f = 4p*+8p'+7p?+4p+1, and L = —72p*+8pt—2p?+9/2. 


The calculation proceeds: 








—4 
— 24/5 














The work to the left of the double line yields 
ut, = 49/2, u, = 1290/8, ug=—3l, u, = 9/2, 
from which we obtain L” from (38), and uj? = —409/12. Thus finally 
J, = —u{/8 = 409/96. 
Of course, the value of u, yields incidentally the result 
J = 49 8. 

In this particular case these results can be easily checked, since on 

inverting #(p) by partial-fraction expansion, we find 
a(t) = fe-#+ ge++e-* cos V7t/4. 

We mention finally an alternative method for both J, and J,. The 
integral for J, in (26) reduces to (3) (with (4) substituted) if g = 0, and 
for non-zero q is of the same form as (3) but with coefficients which are 
polynomials in g. Thus the integral can be evaluated by the method of 
section 2, all elements being, however, polynomials in g. Now since only 
the coefficient of q in the integral is required for J,, only terms of degrees 


0 and | in g need be retained in the elements at every stage. In order to 
carry out the work it is convenient to use the symbol (a,b) for the element 
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(a+6q-+-...). Then we have the rules of operation: 
(a,b)+(c,d) = (a+b,c+4), 
(a, b).(c,0) = c(a,b) = (ca, cb), 
(a,b). (c,d) = (ac,ad+be), 
(a, b)/(e, 0) (a, b)/¢ (a/e, b/e), 
(a,b)/(e,f) = (a/e,b/e—af/e?). 


Similar remarks apply to J,, as given by (34) in which the integral 
reduces to that for J in (11) if g = 0, and is of the same form if q + 0. 
A particular consequence is that the modified-Routh-algorithm method 
for J described in section 3, which cannot be used in conjunction with 
(35) and (36) for J,, since it does not yield u(p), can be used for J, by 
means of the procedure just described. We replace f by (f,f’) and L by 
(L, even part of h’h*), pairing corresponding coefficients together in each. 
Then the algorithm and (22) give, in place of J, the pair (J, —2JJ,). 

We shall illustrate this procedure by using the algorithm to find J, for 
the previous example. 

bel B; 





Yo (9/2, 30) ’ 3, 
fo (1, 4) . (4, 0) (9/2, 12) 
yi > 244 10, 144) 
r; 4,14 (8, 16) (—67/8, —507/16) (1/4, 1/8) 
Ys (57, 1063/2) ” 
r, (5, 19) (4, ©) (57/5, 3149/50) (4/5, —6/25) 
Ys (— 588/5, — 6298/25) 
r; (24/5, 424/25) (— 49/2, 409/12) 


Since f,, is 4 (strictly speaking, the pair (4,0)), (22) gives 
(J, —2J,) = (49/8, —409/48), 


whence we obtain J and J, as before. 


5. Evaluation of /, and J, (r > 1) 


It should be clear that the methods given in section 4 for J, and J, are 
extensible to higher order integrals J, and J,. For example, corresponding 
to (24) and (26) we have 


I, = coeff. of (—q)" in = | #(p)i(q—p) dp (39) 
o 


r 





* = © = 21." /9 4 rlér ’ 
= coeff. of qg” in - h(k+qk' +-q?k”/2+-...+-q" k”/r!) 


a r ” RS (40) 
ari |} figtag' +99" /2+...+4°9"/r!) 
Fs 


+ We use k, g* to denote rth derivatives, in order to avoid confusion with the notation 
u"), v'” for members of sequences of polynomials. 
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so that if polynomials u, u®,..., uw”; v, v®,..., vu” can be found such that 
f(vtqv?+..4-q'v)+(g+¢9' +... +99 /r!)\(ut+quO+...+¢qu”) 

= h(k+qk'+...4-q’k"/r!), (41) 
then =r! ae _slfa- (42) 


Now (41) gives the sequence of equations 


fv+gu=hk=N, (43) 
fv®+gu® = hk’—g'u = N®, (43a) 
fv®+gu® = thk"—g'u—Ig’u = N®, (43 b) 


The first two of these equations have already been met in finding 4,, for 
which it was necessary to solve (43) for the whole of u, use this to determine 
N® in (43a), and finally find w))_, from (43a). Next, to find /,, we must, 
after finding NV“), solve (43a) for the whole of u™, then use wu and u® to 
determine .V® in (43 b), and finally solve this for u‘?_,. This process can 
be continued indefinitely for finding any of the /’s. Once the basic set 
of multipliers and modified coefficients (in the Gauss—Doolittle process) 
has been formed, it is an easy matter to produce J, J,,... in succession in 
this way. 
Similar remarks apply, with appropriate modifications, to the J’s. 
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SUMMARY 


x 


In Part I (1) it was shown that the evaluation of integrals of the forms f xy dt 


@o 
and { tx? dt, where r = 0, 1 


paees and x and y have prescribed rational function 

0 
Laplace transforms, depends on the determination of polynomials u, v sich that 
fv+qu = N, where f, g, and N are prescribed polynomials. A method of solution 


of this equation by means of continued fractions is described which is more economical 
than standard methods. 


1. Introduction 


CONSIDER the integral I | ay dt, 

0 
where x(t) and y(t) have prescribed rational Laplace transforms 7(p) and 
9(p) vanishing at infinity, and such that the poles of Z(p) lie to the left 
of those of 7(—p): write also 


tp) =h(p)if(p), 9(—p) = k(p)/g(p). 

Then, as was shown in (1), there are unique polynomials u and v of degrees 

respectively less than those of f and g, such that 
fv+gu=hk =N, say: (1) 
then I = coeff. of 1/p in u/f. (2) 
In (1) equation (1) was replaced by a set of simultaneous equations in 
the coefficients of u and v, and J obtained by considering the solution of 
this set of equations. An entirely different method of dealing with (1), 
based on the use of continued fractions, will be described in this paper. 


It yields a much more economical method of evaluating 7. A special form 
of the method applicable to the integral 


J f x? dt 
ry 


[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 4, 1959] 
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is also described, as are extensions, analogous to those described in (1), 
for the moment integrals 


I, = | tay dt, J,= | tatdt (r = 1,2,...). 
{ i 


2. Evaluation of J 


Let the polynomials f(p) and g(p) be of degrees m and n, and suppose 
for definiteness that m <n. From the hypothesis concerning poles of z 
and 7, f and g have no common factor. Thus we may express g/f as a 
continued fraction of the formt 

G a 
o isa. see. wl 
f —b,+ p—b,+ 
by means of the h.c.f. algorithm 
fn = GE, + £,, 
a, E, = (p—b spthciiaad 


ay a 
—3- 3 
+ p—b,, @) 


An 1 Ey are. —b »aEn 
Ay Ey = (P—Pn) Ens (4) 
f, the a’s and 6’s are constants, and EL, is of degree 
r = 1, 2,..., so that Z,,,, is a constant. If 
ip) = fothipt--+fmnP™ 
UP) = Jot IPH +9nP"s 


a convenient computational arrangement for the algorithm (4) is the 
following: 


m+. 


here E, : 


g, E, 


m—r-+-1, where 



































row- 
p*p*'p*- we. pp 1 ip 2 op 1p? sums 
Ey g InIn—-1In—2 + 8o 
Ey f SmI mF m—2 +++ 8; 
d IniSn E’ = E,+dp”"E, Gn—r° 
d’ = —9n_1/fm E” = E’+d'p*™"E, th 2° 
qin—m-}) ; Em) ; , . P gh - m) u 
] n—m ~ = 
; E, = E-™+d™)E, ey... 8 
ay €2/fm E, = a, E,—pE, € ++ 
b, —€5/€s E, b E+ Bs Stun 83 
G, = ¢,/€, E; = a,E ~ phy Cp ces 
bs — €3/€s E, b Ey E; Cg oe & 
Om 1 ¢ a em En, > 6... 1 En, + En Cm+t 8mii 
Aan < €m+1 em Dm4t = Om, E,— pE,,, Oma 
b mi mit em4 1 
+ In exceptional cases quadratic or higher denominators may be present in the fraction 
The method is then in general inapplicable. 
5092 .48 Ll 
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The row-sums 8; indicated in the final column provide simple checks on 
the work. If g = d+d'+...+d™-™, 


the following relations, obtained by substituting p = | in (4), must hold, 
and furnish successive checks: 


85 = 8,7 88 


2 0 Bd 


8 = 4,8,+(b,—1)8, 


8 a +-(b,, ,—l)s 


“m? 


—l)s (5) 


m+1° 


m-15m-1 


L(b 


m+1 ~~ 


QoQ = a 


8 
m m m 


Suppose that the successive convergents of (3) are 


Ry OS R, R 


m 
: =F ecug eee 
So I Sy Sin 
Then R,, and S,, are proportional to g and f respectively. Now the R, 
and S; obey the recurrence relations 
R; ; (p -b,)R; 17 a; R; 2 (2 = Zosvey MO); (6) 
S; (p -b,)S; it a,S; : (¢ = 2...., m). (7) 


Since S, = 1 and S, = p—),, it follows from (7) that S; is of degree i and 
its leading term is p'. Thus 

f = SnSm (8) 
and so g = j,i... (9) 
Again, it follows from (4), since Z,,,, is a constant, that a,a,...a,,£, has 
leading term E,,., p”, whence 


Ay Ag... An fy, Ent (10) 
Now, as is well known, (6) and (7) lead to the result 
R,8,;_,—R;_, 8; = (—1)*"a,a,...a; (t = 1,...,m). (11) 


Combining (8), (9), (10), and (11) with « = m then gives 
fRm-1—Bn ee (—1)"E,,41) (12) 
which is the basic preliminary result needed. 


If now u and v satisfy (1), we have by (12) 





NS. 
, sas Sn it uR,, at(—1)°"2,, nz (13) 
from which it follows that if the degree of u does not exceed m—1, 
Uu NS 
MU _ gp|NSm—1) 14 
f ‘| Sn 
where A = (—1)™*!f,, Baas (15) 


and the symbol ‘fr’ denotes ‘fractional part’. 
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We note that (14) always gives a solution of (1), whatever the degrees 
of the polynomials. For if (14) holds, 
NS,, ,/f—(—D""'E,,., u/f = polynomial, say Q, 
and by (12), S,,-11N —ug) = f(Q—uR,,_,). 
But f and S,,_, have no common factor, by (12), so that f divides N—ug, 
i.e. N—ug = vf, say, which gives (1). 

It is clear from (13) that if we impose the condition that the degree of u 
is not to exceed m—1, then (1) has a unique solution, given by (14), what- 
ever the degrees of N and v. In our present application these are at most 
equal to m+n—2 (since # and g vanish at infinity, by hypothesis) and 
n—l. 

By (2) and (14) we see that J is the coefficient of 1/p in the expansion 
of NS,,_,/AS,,, and may be obtained by calculating S,, S,,...,S,, in succession 
by (7), and finding the first term in the remainder when NS,,_, is divided 
by S,,: if this is ap", J = a/A. However, a shorter procedure is available 
which depends on deeper properties of continued fractions than have been 
used so far. These we shall now derive.t 

The essence of the present method is that if the expansion 


SS... c Cc 
mot — 9474... 
P p 


Sin 
is somehow obtained, and if 
hk=N=N4+N,pt+..+N,p” (v<m+n—2), 
then clearly, by (2) and (14), 
I = (Noto+...+N,¢,)/A. (17) 
Thus the essential problem is to determine the coefficients c; in (16) with- 
out determining the numerator S,,_, of the generating fraction (S,, of 


course being known by (8)). 
Now it follows easily from (7) that 


p—b; a; 


—| oa 
—lI 


O 


+ Most of these properties are to be found, essentially, in, for example, Wall (2) where, 
however, both their derivation and their application are presented in a rather more elaborate 
manner. 
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Thus if A denotes the matrix 





“by, —An 0 7 
l o. —An—1 
l 
A : (19) 
—Oe 
a “3 I a 





then S,,_,/S,, is the leading or (1, 1)-element of the inverse of p1— A, where 
1 is the unit matrix of order m: thus, 


S 
, S * = [(p1—A)"]y, (20) 
so that formally we have 
S I l l 
mt | rere At... , (21) 
Sin Pp p~ Pp 11 


It is easy to justify this expansion. For S,, is clearly non-zero and the 
inverse matrix exists for all sufficiently large values of p. Now we have 
identically 


(pt A)(-1 5. eine ie A’) 1—_*_ are, 
Pp p* ow aii 
whence "3 +... Pad = [(p1—A)-"],,+ O(p), 
p z= 11 


since, as is obvious, diagonal elements of (p1—A)-~! are of order p~' and 
other elements are of order p~?. It follows that (21) holds as far as terms 
of order p~’-!, for any value of r. 

Comparing (16) and (21), we obtain 


e, a= {(A’),, (¢ = 0,1,..,). (22) 
Thus if K denotes the m-element column (1, 90,..., 0)’, 
¢ = 1, c, = K’A‘E, 23) 
and if columns K, are calculated in succession by the relations 
K, = K, K, = AK,_,, (24) 


then c, is the leading element of K,. 

In principle this solves the problem stated above. In practice the work 
can be considerably reduced by noting that if y > m, it is only necessary 
to evaluate the c’s in this way as far as c 


m—-1» SINCE Cy, Cyiy--- CAN be 


calculated in succession by means of the recurrence relations 


foCethrCesrt--thmCoim = 9 (8 = 0, 1....). (25) 


ym m 
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These relations may be proved by taking N = fin (1). Clearly the (unique) 
solution is u = 0, v = 1, and (25) then follow from (14) and (16) by con- 
sidering terms in p~*~! in the expansion of 


(fot fi Pt---+Fmn PoP +e/p? + ..-)- 


It is of interest to note that (25) also follow from the Hamilton—Cayley 
theorem applied to A. For since f = f,,S,, = f,,det(p1—A), we have by 


this theorem foit+f,A+...4f,,A™ = 0. 


If we multiply this by A*, and consider leading elements, we obtain (25), 
in view of (22). 

If in (1) we write N = g, so that u = 1, v = 0, we obtain from (14), 
analogously to (25), 


Jor I 91 Oris T eee f I» Crin — (—1)™*1 2, 415) n—1 (r — 0, 1,....m—1), (26) 


where 5 is the Kronecker delta. If vy > n and the c’s are obtained from 
(25), then (26) can be used to check their values. Probably one such 
check, that for r = v—n, would suffice for checking all the c’s. 

The work of calculating the c’s can be reduced still further. In finding 
C1, Cg,-..5 Cp— from (24), only the leading elements of the K, are needed, 
and in view of the special nature of A and K, we may replace A in (24) 
by the submatrix formed by its first [4m] rows and [}(m-+-1)] columns. 
The vectors K, are correspondingly reduced to |}m] rows, and the com- 
plete scheme for calculating the c, up to c,,_, is as follows: 








‘Om, Gm 0 = I Cy Ce . e ° Cmr—1:lm—1 
l b, 1 —a,y, 1 m m “J m Me 
Bie even | odd .~ , r 
‘<: even, 
Nei ; /™ 
See : XL 7 ~odd 
O tie —A4n-+5)) ° Se 
a I ’ bi4(m +301 —%yon+9)1 J (27) 


Here the boundary arrangements in each matrix depend on whether m is 
even or odd. The second matrix has 1’s in its main diagonal and 0’s below 
this. Its columns K, are calculated in succession by (24).¢ Obviously, if 
v < m—1, the scheme will be similar but can be reduced still further, 
using only [}(v+-1)] rows and [4(v+-2)] columns in the first matrix, 
[4(v-+-1)] rows in the second. 

If the scheme (27) is slightly extended so as to yield c,,, the values 


+ Row x column multiplication is used here. If desired, this can be made column x column 
multiplication by transposing A. 
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Cy,.++5 Cp, can be checked by using (25) with s = 0. Further checks, if 
desired, are provided by the relations 


2 2 2 2 
Ci k? 1 —ay, k? Pam Ay 1 k? 5- Am Im -1 4m oki, T 
Coist kis Kiss Am ki k, +1,277 Am Im kis k; +1,3 akg (28) 


where k,, is the sth element of K, (so that c, = k,,).+ It will be seen that 
the scheme (27) contains enough information for applying the checks (28) 
to any of the c, up to c,,_, (and c,, if m is even). 


To illustrate the whole process it will be used to evaluate J for the case 
already treated in (1), section 2, in which 




















. 3+ 2 3n+2 hip 
i(p) Pp * ames i (7 ) 
p'+2p?+3p?+p+1 S(p) 
‘ 4p*— p® + 3p?+3p—1 k(—p) 
ip) = —— : ; 
p? + 3p* + L5p3+- 36p? + 40p4- 27 g(—p) 
so that N (—2, —3, 14, —11, 5, —8, 5, 4, 0). 
The h.c.f. algorithm is as follows: 
p> p® p? p} p® row-sums and checks 
] 3 15 36 40 27 10 
l 2 3 ] l 8 
5 2 37 3f 27 
, 4! 1 f 12 3 39 
eT, est 22 22 —44 22] —22 = 10+(—4).8 
a; bj 
22 66 22 44 22 
3 —88 88 88 88 22.8+(—4).(— 22) 
4 @-—¢es 266 
0 88 8&8 0 = 4.(—22)+(—1).(— 88) 
l 0 88 
0 88 88 1.(— 88)+(—1).0 E 
l 88 
l 








If the c’s up to c, are found from the scheme (27), and ¢,,..., c, (the 
degree of N being 7) by the recurrence relation (25), which now becomes 
= Ls 9 
Cp = —(C,_¢+C,-3+3C,_¢+2,_,), 
we obtain the following scheme (in which elements necessarily equal to 1 

in virtue of their positions have been bracketed): 





Co % Gg Co % Ce Cg Cy 
~— a a ee ee SAT ee oe 
me 
(1) 0 —] ‘ (1) 1 i 
; 7 
(I) 7% 


+ These are particular cases of (53.3) in Wall (2). Incidentally, this can be proved easily 
by means of our fundamental relations (23), which Wall does not give explicitly. 
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The agreement of the values found for c¢, by the two methods forms a 
check on the work up to cy. The remaining c’s may then be checked by 
(26), for, say, r = 2. 

With A -(—88).1 =: 88, we obtain immediately from (17) that 

21/88, as found in (1) by the Gauss—Doolittle method. 

For this example the total numbers of essential entries, multiplications, 
and divisions in the continued-fraction (C.F.) method are respectively (dis- 
regarding checks) about 64, 51, and 15, while in the Gauss—Doolittle (G.D.) 
method they were about 94, 79, and 11. For the general case, leading 
terms in the expressions for these numbers in the two methods are approxi- 
mately: 

C.F. 3m?/4+n?/2, m(n+v—m/2), m+n 

G.D. m?/24-3mn, m(m?+-2n—3m/2), m?/2. 
It is clear therefore that for large m and n the C.F. method is superior in 
every respect. 


3. Evaluation of J 
For J, which is the special case of J with y = 2, so that 
yp) =f(—p) =f*(p), say 
(using the notation introduced in (1), section 4), we must find the coefficient 
of |/p in the expansion of u/f, where 


fut*+f*u = hh* = 2L, say. (29) 


The method of section 2 can be used for this purpose also, but involves 
some wastage of effort since half of the c’s obtained, viz. c,, c3,... are not 
needed, as L is even. We therefore describe now a modification, the 
simplest of several possible, which avoids this waste. 

We split f (which must be a Hurwitz polynomial) into its odd and even 


parts: f=E,+E, (30) 
— Ey = fn P™+Sm-2P™ 2+ 


E, = fins P™ *+Sn-3 P+ 
Then f* =(—1)™4,—£,) = (--1)™(24,—f). 
By means of the algorithm (cf. (4)) 
a, E, = pE,+&,, 
a, E, = pE,+E, 
Amy Egy = PE mt Emsr 
4, Em = PE msi 
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m+ 


where 4,,..., @,, and E,,,, are constants, we obtain the C.F. 


- ee (33) 
m S++ +79 


Suppose as in section 2 that the successive convergents of (33) are 


Then (10) and (11) still hold, while (8), (9), and (12) become 


E, Sa Sa (34) 
E, | = R,, ’ ( 35 ) 
E, R,, a E,S,, 1 ( = 1)" 241: (36) 


The last equation may be rewritten 
fSn—-1— Ey (Sm-1 t+ Rm) = (—D)™ Ena (37) 
while (29) is equivalent to 
fiu*t¥—(—1)™u}+(—1)™2E,u = 20. (38) 
Comparing (37) and (38) we obtain, analogously to (13), 
2L(Sn-1t+ Rina) /f = (Rn +S8,,-1){u* —(— 1) ™u}4 


+(—1)"2uS,,_,+2uE,, .,/f, (39) 
from which it follows that 


u gel Sin it Rin—1)\ 


i 


where o = f,, Bass (41) 


(40) 


- 


- 2m—2, (40) always gives a 
solution of (29), for if uw is given by (40), we have, using (34)—-(36) and the 
fact that the S;and R, are odd or even and of degrees i and i— 1 respectively, 


u u* sos te (4/8 wit R-1 ‘a fr j2L)_ 2h 
all oe \o\ S,+8,, 8,—K,, |) | oo RP ie 
Further, it is easy to show that the difference between any two solutions 
of (29) is expressible as df, where ¢ is an odd polynomial. It follows that 
(40) gives the unique solution of degree less than m (and moreover that 
there is no solution of degree m). 
Now by (34)-(36), 
S__.+R,, , 4... < 1)" Ens a of 2 } (42) 
S,+k S fn Sin Sen +R oo 


m m 


Conversely, if L is even and has degree < 





m) 


t Since f is a Hurwitz polynomial, this form of expansion is always possible, and the a’s 
are necessarily positive. 
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It follows that if S,,_,/S,, (which is odd) has the expansion 


Sin ee fo , © 4% r Com—2 i Com wis 
== — mer T °° - ee —— -  20t*9 
Sin P p* Pp’ p™ 1 pm 
then we may write 


, , 
ma , .. + am=2 , Sam-1 , Com (44) 


1 o™ 437 | pti ae. 
Since L is even and of degree < 2m—2, it is clear that J, the coefficient 
of 1/p in (40), depends only on the coefficients ¢p,..., ¢2,,-2 common to the 
two expansions (43) and (44), and 





where » (< m—1) is the degree of h(p). We shall see later that the expan- 
sion (43) may be used in (40) instead of (44) to calculate not only u 
which yields J, but the whole polynomial u, which yields //,. 
The c’s in (43) are calculated as in section 2. We find 
Co, = K’(A?)'K, 


0 —a, 


l i O 


m~1? 


where now 


l 0 








O 2 eS 


Now the only non-zero elements in A? are in the main diagonal and in 
the parallel lines next but one above and below it. Thus if the vectors 
A°®’K, (A?®)K,... are evaluated they will be found to have zeros in rows 
2, 4, 6,.... It follows that as far as the c., are concerned, A? may be 
replaced by a matrix B obtained from it by deleting rows 2, 4, 6,... and 
also columns 2, 4, 6,..., Viz. 


—a,, Am 4m-1 O 
l —Am-1 —Un-2 Am -24%m-3- 
l —An-3—Am-4 


(45) 
-a 


O : 2 








where a = 4,43, = —a,—a, if m is even, a = a,a,,b6 = —a, if m is odd. 
This matrix has exactly the same form as A in (19), but has [4(m+1)] 
rows and columns. Moreover, to calculate Cg, C4,..., C2 OT C,—, we need 
only use a reduced matrix as in (27), having [}(m-+-1)] rows and [(}m+3)] 
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columns. The remaining c’s, ¢,, OF C,,.1,---, Com-2 May then be calculated 
by means of the recurrence relation 


Smest+Sm-2%s-2++» = 0 (seven > m). (47) 


Of course, if A(p) is of degree less than the maximum, m—1, fewer c’s 
need be calculated. 

The h.c.f. algorithm is similar to that in section 2, but simpler, having 
only the lines for Z,, E,,..., E,,,,. It may be checked by row-sums in the 
same way. To check the c’s up to c,,_, or ¢,,_;, We may use B to calculate 
C,, OF C»,,, and check it by the recurrence relation (47). Alternatively, we 
may use the following relations, which replace (28) in this case: 


. 2 2 3 
C4; kei T Am Am ASi2 T Aan Im 14m 24m 333° aeey 
Case = Kei heise, t+ On Im—thei eiseget+--- (48) 


Finally, an overall check, replacing (26), is given by the relations 
: g ; 


fm—1Cst+JSm—3Cs-2t--» = (—1)"*"E,, 1 5¢0m-2 (m—1 < s(even) < 2m—2). 
(49) 
To illustrate, we shall evaluate as in (1), section 3, J | x* dt, given that 
0 

=(p) p*— p®+ 3p*— 2p?—5p?+p+3 h 

x(p 7 . A ae . , 

F p' + p*+ 6p®+ 5pt+-9p3+-6p?4+-3p+i ff 

so that L = }(9—31p?+-47p*— 26p*— 5p§+ 5p!+-p"?). 


The h.e.f. algorithm is as follows: 




















row-sums and 
a, Pp’ p® p> p* p® p? p} p° checks 
| 4 6 9 3 19 
l 5 6 ] 13 
l 
l 3 2 6 1.19—13 
l 
2 4 ] 7 1.13—6 
2 
2 3 5 2.6—7 
1 
1 1 2 1.7—5 
$ 
p= 4-5-2 
4 ” 
, E,= 4} $= 3.2-4 
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We now find ¢,, ¢y, ¢g, (¢g) from the reduced matrix B: 
& 2 | qe GQ “w “ts 
\@) -—1° 9 = ‘—9 21 


4 


A 


8 
TA 





and ¢,, Cy», Cy, from the relations 

Co = — (6c, 2+ 9, at3ey,_¢), 
an overall check being provided by (49) for, say, s = 12. Then by (41) 
and (45) we have o = 4, and J = 129, as found in (1). 

In the general case, with A(p) of maximum degree (m—1), the whole 
process involves roughly 5m?/16+3m entries, 5m?/8—m multiplications, 
and 3m/2 divisions, and, by comparing these numbers with those given 
at the end of (1), section 3, is seen to be superior to either of the methods 
described there, except perhaps for m > 24, when the number of multi- 
plications and divisions is greater than in the method based on Routh’s 
algorithm, though the number of entries is far smaller. 


4. Evaluation of /, and J,,r > 1 
It was seen in (1), section 4, that 
Ty = U.S (50) 
where u‘!), is the leading coefficient in u, determined from the pair of 


m—1 
equations futgu = hk = N, (51) 
fv&+gu® = hk’—g'u = N®, (52) 
Thus to find J, it is necessary to solve (51), which is identical with (1), for 
the whole polynomial 
U = Up+uU, p+...tu,,,p"", 
and not merely its leading coefficient u,,_, as suffices for 7. Then, having 
substituted u into the right-hand side of (52), which is of the same form 
as (1), we obtain J, from this in the same way as 7 was obtained from (1). 
Now u/f is given explicitly by (14), and we shall now show how the c’s in 
the expansion (16) can be used to determine u in the simplest possible way. 
We have, by (14) and (16), that 


yiCy , © c c Cy 
= Sot hiv t hav) N(2+ S54) +0, (24 3 4.) 4-0,(2 +f], 


P 
(53) 
If we use the fact that for all sets of values of N,, N,,... the right-hand side 
of (53) must be a polynomial, i.e. its fractional part must be zero, we 
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immediately recover the recurrence relations (25). It will be convenient 
to use the notation 


bin = elo tSerr Ceri t thn Orsm-i: (54) 
The equations (25) may then be written 
top = 0 (r= 0,1....). (55) 
We note also that t= f_¢,. (56) 
Considering now terms in p*~! in (53), we obtain 
Au;_, = s tN, (8 = 1,3.,..., 9). (57) 
r=0 


Now (54) may be written in either of the two forms 
(a) bssarst tip fie, 
(b) ti, t; Hirsi fie, (58) 


These imply that the matrix T of elements ¢;, (¢ = 0, 1,...,m,r = 0, 1,...) 
may be built up by using fo, f;,..., f 


m 


to label the rows, Cy, ¢,,... to label 
the columns, and, starting with known element values, moving (a) 
diagonally down and right, or (b) diagonally up and left. In a step (a) we 
subtract the f—c product corresponding to the initial position; in a step 
(b) we add the product corresponding to the final position. 

For steps (a) the starting values are foo, ty;,... Which by (55) are all zero. 
Initially, only Cp, ¢,,..., C»—, are assumed known. By a succession of steps 


(a), starting with fy, and using Cp, ¢),..., C,—,; IN succession we arrive at 
t 


ty, and using ¢,, Cy,..., ¢ 


mm Which in turn by (56) yields the next of the c’s, c,,. Starting now with 


_ we obtain f,, .,4, and hence c,,,,. Continuing this 


” mm 
process we can obtain all the ¢;, with r > @.t 
To obtain the remaining t’s, we start with t,,9, tinys> tmom—1 given by 


(56), and use steps (6), i.e. move upwards. The arrangement of the work 
is indicated in the following table: 





Co Cy Ce Cm—1 Cm C, 
fo 0 0 0 0 0 0 0 
4 4 4 4 “ a 
fi tio ty tie th, 
. 
. . 
(59) 
Im 2 tm 2.0 
Im 1 tm 1,0 tm 1,1 
* * * 
4 ~ 
Su tmo tiny tm $ . ’ tmm 1 tam ° . . , bm 











+ Note that the procedure given in § 2 for finding ¢,,, Cm 4,5-++ 


.e. using (25), amounts to 
performing the subtractions involved in these downward steps, and so obtaining typ, 
tm.m+is-++» tmy, Without however recording the intermediate ?’s. 
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Two relations connecting successive t’s, which may be used to check 
the completed table, can be obtained by giving the N’s special values. 
Taking VN = p*f, so that N, = f,_, (r > k); 0 (r < k), we obtain from (1), 
as the unique solution with u of degree less than m, u = 0, v = p*, and 
so by (57), 


m 


Dbizwle =O (k =0,1,..., § = 1,2,...,m). (60) 
s=0 


Thus the full table (59) may be checked by applying (60) to successive sets 
of (m+-1) consecutive t’s in each row. 
By taking N = p*g (k = 0,1,...,.m—1), we obtain similarly u = p*, 


v 0, and 
n 


> tiirGr — VB; eat (k = 0,1 grees m—1). (61) 


r=0 
(61) does not hold for k > m—1, when u is of degree > m. 

Having found the t’s, the u’s, multiplied by A, are found from (57) by 
multiplying the various rows of the table by the vector N. The right-hand 
side N“ of (52) can then be found. This has degree m+n—2, and if 
v < m+n—2, it will still be necessary to find ¢,,,,..., C,.,—2 by using the 
relations (25). Then J, is given by 


I, = (N@eg +N, +... +N, sman—s)/A- (62) 


To illustrate the process, let and 7 be the same as in the example of 
section 2, where m = 4, n = 5, vy = 7. Using the values of ¢y, ¢,, Cg, Cs 
found earlier, our table becomes: 











want 1 o —-1 O 2 -3 1 u, 
fr ™ 

; Tt 2 0 0 0 0 0 0 0 

l 3 -1 —1 0 1 0 —2 3 -10 
3 5 2 2 —] 1 1 —2 1 — 42 
2 3 2 l -2 2 1 —§ 7 1 
¥ l 1 0 Ls Ss 2. -3 1} -21 
N, 2 3 14 1] 5 8 5 4 











Since A 88, we have 


u = (—10—42p+ p?—21p)/88, 


and we find 
N® = (—928, 888, 1294, —394, —1421, 505, 1665, — 105, 0)/88. 
Since m+n—2 = 7 = v, there are no additional c’s to be found in this 
case, and we have finally, by (62), J, = —467/968. 
Although the extra work involved in finding J,, as compared with J, is 


greater than in the G.D. method, the total amount of work remains 
smaller. 
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Turning now to J,, we have, as seen in 1, (37), 


J, = —ul_,/2fras (63) 
u\)) , being determined from the pair of equations 
fu®*+ftu = hh* = 2L, (64) 
fur* tf *y) . 2L, (65) 
where L” = even part of (h’h*—f'u*). (66) 


Thus when L” has been found, J, is given by 
J, (LM eg+ LW egt... + LD _ 9 Com—2)/20f,n- (67) 
To find L® we must solve (64) for u, which is given by (40). As at (57) 
above, we may write : 
a 
ou,.,= >t, L, (¢ = 1,2....,), (68) 
r=0 
where 2, the degree of L, is < 2m—2, and the ¢;, are again given by (54), 
though it is to be noted that the c’s are the coefficients in the expansion 
(44), not (43), and may be needed beyond ¢,,,_, (if 2) > m—1), whereas 
(43) gives them correctly only up to ¢,,,-5. However, this is immaterial 
since only the values up to c,,_, are needed for starting the table (59), 
which is then completed by upward and downward steps as for J. 

A simplified version of the table can be used for the present purpose. 
In the first place, since L is even, it is only necessary to find the ¢;, for 
even r up to 2;t moreover, since c, = 0 for odd r up to 2m—3, 

tere = bisarar (FP 1, 3,..., 2m—3, t = 0, 1...., m—1l), (69) 


so that in the complete table of the ¢;, (r < 2y < 2m—2), nearly half the 
entries are repetitions. Thus an abbreviated table may be used as follows: 








€o C2 C4 
ho 0 0 0 
hi tio 0 0 a 
» 
4 4 
fe too tos tog a 
4 
fs ts tye tos ry . . . = 
A \ ( #0) 
4 \, 
a 4 
Sun 2 bmn 2,0 ’ ~ 7 . * 
K, * 
te 1 tm1.0 tm 1,2 
Su tno tn 2 tm 4 








+ As before, if 2n < 2m —2, it will still be necessary to find additional c’s, up to Cy,» 
for insertion in (67). 
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The procedure is the same as before except that the steps are one place 
to the right (or left) and two places down (or up), and we have not only 
too = tog = «.. = 0, but also t,, = t,, = ... = 0, these providing starting 
values for all the downward steps. Moreover, on the two bottom lines of 
the table, to start the upward steps, we have 


true =Sm-i1% (r = 0,2,....m—2 or m—3), 
as well as tar =Smc, (¢ = 0, 2,....m—2 or m—1). 
To check the table, we use the relation 
tm—-ue =Sm-1%r (7 = m or m—1.,..., 2m—6, 2m—4). 
Relations analogous to (60) and (61) can be obtained by taking 
(i) 2L = p*(f+f*), keven, << m: u= p* = u*, 
(ii) 2L = p*(f—f*), kodd, <m: u= —p*, u* = p*. 


Using (68), these give 


bifottirrafet--- = 08:54, (6 = 1,2,...,.m,r = 0, 2,...,(< m)), (74) 


tirlittirsefst-.. = —od, (t = 1, 2,...,m, r = 2,4,...,(< m)), (75) 


which may be used as additional checks. (Equations (73) do not suffice to 
check the whole table). 
To illustrate the process, we take the example in (1), section 4, in which 


E(p) = (12p+- 16p?+ 10p-+-3)/(p+1)(2p+-1)(2p?+- p+1), 
and =f = 4p*+8p3+7p?+4p+ 1, L = —72p*+ 8p*—2p?+3. 


The calculation proceeds thus: 








4 7 1 Co Cy Cy C6 
1 —1/ 1/10 —1/8 


0 0 0 0 
12/5 0 0 0 
31/6 —1 / —1/10 

8 —8/5 / — 2/5 

4 —4/ / —1/2 


9/2 —2 ‘> ae 
12 —233/4 49 





























Cy Cs 
3/50) \Q) —1/5.°1/107 
4 - 


‘ F 
Sed 
\ 


from which cu{})_, = —409/5, and J, = 409/96, as before. 
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The evaluation of J, and J, (r > 2), using the C.F. method, involves an 
extension of the processes above, indicated by (42) and (43) of (1), and 
need not be further dealt « ‘th here. 
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